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ABSTRACT 


Modeling studies of the interaction of a delta-wing 
aircraft with direct lightning strikes have been carried out 
using an approximate scale model of an F-106B. The model, 
which is three feet in length, is subjected to direct injec- 
tion of fast current pulses supplied by wires, which simulate 
the lightning channel and are attached at various locations on 
the model. Measurements are made of the resulting transient 
electromagnetic fields using time-derivative sensors. The 
sensor outputs are sampled and digitized by computer. The 
noise level is reduced by averaging the sensor output from ten 
input pulses at each sample time. Computer analysis of the 
measured fields includes Fourier transformation and the compu- 
tation of transfer functions for the model. Prony analysis is 
also used to determine the natural frequencies of the model. 

Initially, electric and magnetic field measurements were 
made with the F-106B model replaced by a simple circular cylinder. 
Analysis of the waveforms has produced important new information 
about the effect of wire attachments on scatterers. The 
natural frequencies of cylinders with wires show much greater 
damping and slightly higher resonant frequencies than do those 
of isolated cylinders. The former is duo to the wires carrying 
energy from the cylinder, adding to the radiation damping al- 
ready present. The latter is due to the wires changing the 
geometry of the ends, which reduces their capacitance and makes 
the cylinder electrically shorter. The results concerning the 
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natural frequencies of the cylinders also give confidence that 
the basic measurement and analysis technique is correct. 

Comparisons of the model transfer functions with 
spectral amplitiMes computed for in-flight data obtained on 
the F-106B by N.A.S.A. show good agreement regarding the 
frequencies of the resonances. The frequencies vary somewhat 
depending on the wire attachment points. The average values 
are 7.5, 13, 19, 24, 29, 35, and 41 MHz. These frequencies 
are all matched in in-flight data except for 19 MHz. In flight 
it is 21 MHz. 

Comparisons of model natural frequencies extracted by Prony 
analysis with those for in-flight direct strike data usually 
show lower damping in the in-flight case. This is .indicative of 
either a lightning channel with a higher impedance than the 
wires on the model, only one attachment point, or short streamers 
instead of a long channel. There is also some in-flight data 
for nearby strikes. Here the damping should be and is lighter 
than on the model. For the first natural frequency (7.5 MHz) 
the normalized damping rate is about -0.16 for the nearby strikes 
and -0.24 for the model. 



chapt:sh X 
INTRODUCTION 

The work presented in this report was carried out as 

part of the N.A.S.A. Storm Hazards research program. It 

represents one step toward a general understanding of the 

electromagnetic environments encountered on aircraft during 

lightning strikes. The work involves ix controlled laboratory 

* 

experiment using an aircraft model and a comparison of the 
resulting fields with those observed on the N.A.S.A. F-106B 
research aircraft during lightning strikes. Exterior electric 
and magnf^t<’c fields have been studied. 

The central issue of the report is the nature of the 
resonances, or natural frequencies, of the aircraft with an 
attached lightning channel. Resonance is a fundamental physical 
phenomenon and represents one of several aspects of an aircraft- 

f 

lightning interaction which are important areas of study. In 
the in-flight situation, the resonances of the aircraft modify 
the fields imposed by the lightning by enhancing the spectral 
components lying at or near the resonance frequencies. Thus, 
some knowledge of the resonances is necessary for interpreting 
the measured lightning data. In addition, the lightning channel 
plays a role in determining the characteristics of the resonances 
including their frequencies and damping rates. Thus, detailed 
studies of the resonances observed during lightning strikes can 
yield information about the channels. The resonances are not 
the same, for example, as they would be for the case of a 
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nearby lightning flash or a nuclear EMP, where there is no 
oh/:;iinel attachment. Of course/ the channels have time- 
dependent and nonlinear properties; our technique does not 
attempt to model these. 

The laboratory measurements were done in the time domain 
using short pulses. The pulses are applied on a wire/ which 
represents the channel . The technique is similar to that dic- 
cussed in a Commission A (Electromagnetic Metrology) session 
at the January, 1982, URSI meeting [1] » It appears that ours, 
is the only time-domain facility where direct current injection 
is used instead of an incident wave input. Some results from 
other workers are contained in the references [2,3,4]. 

Our approach to the evaluation of resonances employs 
two types of analysis: first, ordinary Fourier spectrum 

analysis and, second, Prony analysis. These and other tech- 
niques have been compared in a recent review paper [5] . The 
application of the Prony analysis [6-11] derives from the 
singularity expansion method (SEM) of Baum [12] where the fields 
on an object are expressed in terms of a set of natural fre- 
quencies. Prony 's method extracts these natural frequencies 
from the time-domain waveforms. The concept of natural fre- 
quencies has already been used to characterize the response 
of an aircraft to an EMP [27] and in radar target identification 
schemes [28,29]. 

The present investigation is an extension of work reported 
previously which consisted of the first results from the modeling 
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experiment [13] and preliminary analysis of in-flight data 
from the P-106B [14] , A few of the results from this pre- 
vious work are included in the present report/ iae the sake 
of completeness. Several related reports, which are not 
referenced in the previous work, have also appeared [15-18] . 

The new contributions of the present work include labora- 
tory measurements of the fields on circular cylinders and on 
the aircraft model for 4 different wire attachment configurations 
Also included is extensive use of Prony analysis on the wave- 
forms from the cylinders and the aircraft model and also from 
the F-106B for both the 1980 and 1981 flights. The presentation 
is as follows: Chapter II describes the laboratory measurement 

techniques. Chapter III describes the Prony (and Fourier) 
analysis, which is applied in Chapters IV, V, and VI. Chapter 
IV covers the cylinder results, and Chapter V gives the results 
for the aircraft model. Chapter VI compares the model results 
with those from the F-106B. Conclusions are drawn in Chapter 
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CHAPTER II 


EXPERIMENTAL AND MEASUREMENT TECHNIQUES 

Apparatus 

The modeling apparatus is shown in Figure 2.1. It 
consists of a ground plane (12 by 12 feet) , p'ulse generator 
(Tektronix 109), sampling oscilloscope (Tektronix 561A*) , and 
a computer for digitizing and recording the data (DEC PDP 
11/04) . The object under test is located 10 feet above the 
ground plane, with wires attached to simulate the lightning 
channel. An approximate scale model of the F-106B, equipped 
with small B-dot and D-dot sensors for measuring the transient 
electromagnetic fields, is shown in Figure 2.1 [13]. 

A pulse of 0.75 ns duration and risetime of 0.25 ns 
is launched at the bottom of the lower wire. The pulse 
travels up the lower wire, over the model, and on up the 
outer conductor of the upper wire, which is actually 0.141 
semirigid coaxial cable. The sensor outputs are carried 
on the inside pf this cable. No measurable leakage to the 
inside of the cable has been detected. This is as it should 
be since the skin depth for the pulse is two orders of magni- 
tude smaller than the cable wall thickness. When the top of 
the apparatus is reached, -a transition is made to 0,5 inch 
diameter cable to complete the run back to the sampling 
oscilloscope, which has an equivalent risetime of 25 ps. 


*The following plugins are used: Type 3S2 with S-4 sampling 

heads and Type 3T2 . 
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APPARATUS FOR 
AIRCRAFT-LIGHTNING MODELING 



Figure 2.1. Modeling apparatus 
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The clear timOj or the time before any reflections return 
to corrupt the data, is 20 ns. This corresponds to the 
time required for the pulse to reflect from the bottom of 
the model, travel back down the lower wire and back up to 
the model after reflecting from the ground plane. 

Data Acquisition System 

A detailed block diagram of the data acquisition system 
is shown in Figure 2.2, The components located between the 
sampling oscilloscope and computer system were omitted from 
Figure 2.1 for simplicity. The computer system consists of 
the computer and floppy disk drive, an A/D converter, a pro- 
grammable clock, a video alphanumeric and graphics terminal 
with hard copy unit, and a printing terminal. The A/D con- 
verter and programmable clock are special boards which plug 
directly into the computer backplane. They are both under 
software control using Fortran subroutines [19] . The A/D con- 
verter is a Data Translation DT 1712 which has one 12 bit A/D 
‘converter with 8 differential input channels multiplexed into 
it. The maximum throughput rate is 35 kHz. The input voltage 
range is +^10 volts, so one least significant bit is 4.88 mV. 
The A/D board requires an external trigger or clock signal to 
mark the beginning of each conversion. This is supplied by a 
DEC KWll-K programmable clock. This board contains a 1 MHz 
crystal controlled oscillator with programmable divider cir- 
cuitry to produce different clock rates. 

The clock can be started by an external trigger pulse, 
which is the basis for control of the data acquisition 
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Figure 2.2. Block diagram of data acquisition system 
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aystem. The oscilloscopo sweep sawtooth is usefl to * mer- 
ate the "data window" with the trigger and window generation 
circuitry shown in Figure 2.3. The leading edge of the 
window is used to trigger the clock. The data window thus 
determines when data acquisition on a particular sweep will 
start and how long it will continue. The window actually 
lasts slightly longer than the 20 ns clear time to make certain 
that the required number of samples has been taken. The window 
is also used to measure the actual sweep time of the sampling 
oscilloscope using the 1 KHz crystal oscillator and counter, 
and to measure the sweep length by observing the window on the 
monitor oscilloscope. Since the sampling oscilloscope assembles 
a waveform from many consecutive signals, its display is not in 
real time and the actual sweep rate in s/div is needed in order 
to compute what the clock rate should be to produce the desired 
equivalent sampling rate. The monitor scope also displays the 
waveform to be sampled and visually insures that it is entirely 
within the sampling window. 

The oscilloscope displays the sensor signal at an equiva- 
lent rate of 2 ns/div and has 10 horizontal divisions cover- 
ing the 20 ns clear time of the experimental setup. This is 
to be sampled at an equivalent sampling rate of 20 GHz, cor- 
responding to a sampling interval of SO ps. This high samp- 
ling rate will remove any possibility of aliasing in the data. 
Since the oscilloscope is sweeping at an equivalent rate of 
2 ns/div, 40 samples/div are required for a sampling interval 
of 50 ps. If the actual sweep rate of the oscilloscope is 
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Figure 2.3. Trigger and window generation circuitry. 
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measured to be T sec/div, then the sampling rate of the 
A/D converter should be 

R =3 40/T samples/sec. 

In order to achieve a jitter free display# the sampling 
oscilloscope itself is set internally to run at the fixed 
rate of 100 samples/div. This gives an actual sweep rate 
for the oscilloscope of approximately 0.13 sec/div and an 
actual sampling rate of about 308 samples/sec. 

The vertical outputs of the sampling oscilloscope have 
a 10 k ohm output impedance which is too high to allow proper 
settling time for the input multiplexers. The settling time 
is the time for the charge injected into the input lines by 
the solid state multiplexers to settle to less than 1/2 of 
one LSB of error. This is a function of the cable capacitance 
and the source impedance time constant. So, if a high source 
impedance is present, only very short input cables can be used. 

To alleviate this problem, buffer amplifiers of the type shown 
in Figure 2.4 have been designed and constructed. They are 
double ended on both input and output. They have a high input 
impedance to avoid loading the oscilloscope and a low output 
impedance to solve the settling time problem. They also are used 
to adjust the gain of the signal from the oscilloscope in order 
to utilize the entire +10 volt dynamic range of the A/D converter. 

Data Acquisition Code 

The data acquisition code written in Fortran is given in 
Appendix I. The code actually acquires 11 successive sv/eeps 
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from the oscilloscope. The first is ignored to allow the 
triggering of the sampling sweep to stabilize. The last 
ten are averaged point by point in order to reduce the base- 
line noise in the measurement. This averaged waveform is 
then stored on magnetic disk. Detailed operator instructions 
for running the code are given elsewhere [20] . 

Model and Sensors 

The approximate scale model of the P-106B delta-wing 
aircraft is constructed as follows; The fuselage consists 
of a 3 foot length of aluminum cylinder, 4 inches in dia- 
m.eter. Flat end caps are machined to fit into each end and 
are secured with screws. The wings and tail are cut from 
1/16 inch brass sheet to scale with the aluminum cylinder. 

The overall scale of the model is 1/18.8 of the actual F-106B. 

A comparison of the model to the actual aircraft is shown in 
Figure 2.5. The agreement is quite good except in the cockpit 
area. No attempt has been made to include the cockpit on the 
scale model and the effect of its omission on the external 
response of the model is expected to be slight. 

Two B-dot sensors and one D-dot sensor have been placed 
on the model in the locations corresponding to their actual 
locations on the F-106B. The B-dot sensors, one sensitive 
to longitudinal currents and one to transverse currents’*', are 
on the fuselage over the starboard and port wings, respectively. 


We follow the terminology used by N.A.S.A. for the sensors on 
the airplane: "longitudinal" and "transverse" refer to the 

direction of current, not the direction of B. 
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The D-dot sensor is under the nose or chin of the aircraft. 

The output polarities for these sensors are shown in Figure 
2.6. Current increasing fore to aft produces a negative re- 
sponse from the longitudinal B-dot sensor, as does current 
flow from starboard to port for the transverse B-dot sensor. 
Electric field increasing in the outward normal direction 
produces a positive response from the D-dot sensor. A photograph 
of the model in position above the ground plane is shown in Fig- 
ure 2.7. Wires are attached to the forward and rear fuselage 
as was shown in Figure 2.1. The two B-dot sensors can also be 
seen in the photograph. 

The B-dot sensor is shown in Figure 2.8, along v:ith its 
equivalent circuit and its response to a step magnetic field. 

It is a loop made by bending 0.141 inch diameter semirigid 
coaxial cable into a 0.9 cm mean radius semicircle and cutting 
a gap in the outer conductor. The size of the sensor has to 
be chosen large enough to obtain adequate sensitivity but small 
enough so the bandwidth will be adequate. It should also be 
small with respect to dimensions on the model. This radius, 

0.9 cm, is the smallest radius in which the 0.141 cable can 
be bent without damage. Simple current estimates for the model 
indicated that this should be large enough for adequate signal. 
The bandwidth, f^, is calculated from the inductance by assuming 
the simple equivalent circuit in Figure 2* 8b, loaded by the 50 
ohm line. The output falls by 3 dB when 2 it f^ L equals 50 ohms, 
which is the expression for f^. Also, the sensor risetime can be 
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Figure II-7. Photograph of F-106B model above the ground plane. 
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For B- Bo u(t) 

V8(t) = ABBo 3(t) 


a. Sensor b. Equivalent circuit 



L 


c. Output voltage for step d. Integrated output 
input voltage 


Figure 2.8. B~dot sensor and associated circuit and 
waveforms. 
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calculated using the formula t^^f^ = 0.35. The inductance of 
the semicircular loop can be calculated by using the formula 
[ 211 : 

L S I 2l| (in 7, 

2 10^ r 4 

where 

g = mean diameter of loop (cm) 
r = radius of wire (cm) . 

The factor 1/2 is put in since the sensor is a semicircu- 
lar loop on a ground plane. Here, with r = 0.1791 cm 
(0.0705 in) and g = 1.8 cm (0.7 in), L = l.l’lO"® H. Thus, 
the bandwidth is f^ = 50/2 ttL = 724 MHz, and the risetime 

is t = 0.35/f = 484 ps. 

r o 

The sensitivity (or equivalent area) of the B-dot sensor, 

A_, can be approximated by the geometrical area of the loop. 

D 

2 -4 2 

The area is 1/2 (tt r ), where r = 0.9 cm, so = 1.27*10 m . 

The D-dot sensor is shown in Figure 2.9, along with 

its equivalent circuit and response to a step electric field. 

It is simply a 1 or 2 cm monopole antenna with a small ground 

plane. It was determined that these lengths would provide 

adequate sensitivity using some electric field estimates on 

the model. The raonopole needs to be short enough, though, so 

that errors due to the fact that the field falls off as 1/r 

from the cylindrical surface will be negligible. The same is 

true of the B-dot sensor described above. Indeed, these errors 

are slight as determined by calculations of correction factors 

due to the 1/r field variation. Unlike the B-dot sensor, neither 

the effective area nor the bandwidth can be calculated for this 
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D-ilQt sonsejr* Tho offoqtivo length can bo assumed t© bo equal 
to tho physioa.1 length, but this is related to the offoetive 
area by tho oapaeitanoe whieh eannot be ealeulatod. 

In order to oalibrato tho sensors, measuring tho effoo- 
tivo areas and risotimos of eaeh, tho sensors may bo oxpoaod 
to a B- or D- field with a stop change in time. The sensor 
output resulting from the stop is integrated and or 
found from the final value roached by this integral. This 
is shown in Figures 2,8 and 2.9, where and aro the 
values reached by the magnotio and electric fields at tho 
sensors. Tho equations aro almost identical in the two 
cases. In Figure 2.8, the B-fiold is given in terms of the 
unit stop, u(t), and the sensor voltage source then contains 
an impulse function, vMt) , Tho sensor output voltage computed 
from the equivalent circuit is shown in Figure 2,8c. The in\- 
pulse charges tho inductor instantaneously, followed by an 
exponential decay. The final value of the integral, from 
Figure 2.8d, is Ag from which is obtained. The same 
procedure is applied in a dual manner to the D-dot sensor in 
Figure 2.9, o.xcept that the final value of the integral is 
Ap R, The sensitivity of the D-dot sensor is then A^ R 
v/(A/m^) . 


The bandwidth of each sensor can be computed by meas- 
uring the 10% to 90% risetime of the integrated output, and 
then applying the formula t f » o.3S. 
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The sensor calibration is done on the conical trans- 
mission line calibrator shown in Figure 2.10. It consists 
of a conical section of transmission line formed by the 
conical monopole and the 12 x 12 foot ground plane. The 
electromagnetic fields produced on this line can be calcu- 
lated by using the formulas [221; 

UqI EoTiV 

* 2iTrsin0 D0(r,0) = iiirKsine 

where 

V = voltage on cone 
I = current on cone 

r = radial distance to point of interest 
9 - polar angle to point of interest 
n = impedance of free space = 377 ohms 
K = impedance of conical line - ri/2*ir Jin cot i{j/ 2 
s cone half “angle. 

The sensor under test is placed on the ground plane in the 
location shown, which gives maximum clear time, 8 ns. A 
step is applied to the cone producing a spherical electro- 
magnetic wave with fields described by the above equations. 
The step is shown in Figure 2.11. It has a very fast rise 
to approximately 90% and then a slow increase after that 
which seems characteristic of line type pulse generators. 

The output and integrated output for the B-dot sensor 
are shown in Figure 2.12. For this measurement, the sensor 
was placed only 2 feet from the cone rather than 4 in order 
to obtain adequate signal. The sensitivity of the sensor 
is determined by the final value of the integrated output, 
as described above. This value is equated to Ag B^, where 




TIME (NS) 


Figure 2.11. Voltage step applied to conical 
calibrator . 




TIME 


Figure 2.12. Output and integrate 
B-dot sensor. 



output for the 
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Bq is the value reached by the magnetic field at the sensor. 

. . -7 

For the positron 2 feet from the cone, - 1.97 *10 T. 

-12 

The final value of the integral is 21.9*10 vs, giving 

-4 2 

the effective area = 1.11*10 m . This is somewhat less 

Jd 

than that computed using the geometrical area. This indicates 
that for this type of B-dot sensor, the effective area is about 
14% below the mean geometrical area. 

In order to account for the risetime of the pulse generator 
and measurement system, the risetimrs of the sensor can be esti- 
mated by using the following method [23] ; 

4 -4 * 4 

^total system sensor 


so 


= / 


/' 


- t 


sensor 


total 


system 


From Figure 2.11, the system risetime is 350 ps. From 
Figure 2.12, the total risetime of the B-dot sensor is 670 
ps , giving an sensor risetime of 570 ps and a bandwidth of 
f^ = 615 MHz, lower than that predicted by the inductance 
calculation above. 

The outputs and integrated outputs for the 1 cm and 2 
cm D-dot sensors are shown in Figures 2.13 and 2.14, respec- 
tively. Similar calculations to those above for the B-dot 

-3 2 

sensor yield a sensitivity of 8.89*10 V/ (A/m ) for the 

-3 2 

1 cm sensor and 22.8*10 V/ (A/m ) for the 2 cm sensor. As 
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noted earlier, the effective length of the monopole is re- 
lated to the effective area by the capacitance as follows: 


If the effective length is assumed to be the physical length, 

then the capacitance can be calculated using the above formula 

and the effective area calculated above. If this is done, 

-13 

the values are C = 1.54*10 F for the 1 cm sensor and C = 

-13 

2.1*10 F for the 2 cm sensor. The bandwidth of the sensors 

can then be estimated in a manner similar to that used for the 

B-dot sensor. The bandwidth is calculated by assuming the 

equivalent circuit shown in Figure 2.9b, loaded by the SO ohm 

line. The output falls by 3 dB when l/(27rf^C) equals 50 ohms. 

The risetime can then be calculated using t f =0.35. This 

r o 

yield'-j a risetime of 16.9 ps for the 1 cm sensor and 23.1 ps 
for the 2 cm sensor. These are more than an order of magnitude 
below the system risetime so they cannot be measured on the 
calibration system. Thus, the bandwidth of the model measurements 
will not be limited by the D-dot sensors but only by the B-dot 
sensor. 

Time Domain Reflectometry 

The technique of Time Domain Reflectometry (TDR) can 
be used to learn a considerable amount about the response 
of the F-106B to an electromagnetic transient. The tech- 
nique consists of launching a fast risetime voltage step 
along the input line to the model and observing the reflec- 



tions with a sampling oscilloscope. The TDR display is 
then a plot of the magnitude of the reflection coefficient 
(p) versus time. The equipment used here is the Tektronix 
7603 oscilloscope with the 7S12 TDR/sampler plug-in. The 
7S12 uses the S-52 pulse generator head (t * 25 ps) and 

Jm 

the S-6 sampling head (t^^ » 30 ps) * The results from this 
technique will be presented in a later chapter. 


CHAPTER IX I 


NUMERICAL TECHNIQUES 

The numerical methods used in analysis of the data 
measured in the laboratory and the in-flight data are the 
fast Fourier transform (FFT) and Prony analysis. A stick 
modeling technique for the aircraft is also used in an 
attempt to calculate the natural resonances of the F-106B 
for comparison with measurements. 

Fast Fourier Transform 

The FFT performs a fast calculation of the discrete 

* M W Gk Ma W W A^ MA k • W 

^o' ^1' * * * ' ^N-1 following summation: 

P„ = ^ f ^-j2TinK/N foj. X = 0,1,..., N-1 

^ n=0 " 

Similarly, the inverse DFT is calculated using: 


f = I F^ gj2TTnK/N ^ _ 0,1,..., N-1 

" KiO * 


The finite data sequence to which the FFT is to be applied 
may be viewed as being obtained by windowing an infinite 
length sequence with a cutoff function. Multiplication 
in the time-domain by this window function means the re- 
sulting transform is the convolution of the desired trans- 
form with the transform of the window function. In the 
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case of a rectangular window, sidelobes characteristic of 
the (sin irf)/7rf function are introduced. Therefore, the 
implicit assumption that the data outside the window is 
zero introduces leakage or power spreading into adjacent 
frequency regions. 

The data windowing also determines the frequency reso- 
lution of the FPT. The resolution is limited to the main 
lobe width of the window transform. For the rectangular 
window, the width of the main lobe of its transform, (sin 
TTf)/TTf, is approximately the inverse of the observation time 
in seconds, or 1/KAt Hz. 

Zero padding, or adding zeros to the end of a data 
sequence before transforming, will serve to interpolate the 
spectral values between those which would have been obtained 
otherwise. In this way, the appearance of the spectrum can 
be smoothed. Also, padding can help to resolve ambiguities 
in the spectrum and to reduce quantization error in the 
estimate of frequencies of spectral peaks [5] . 

The FFT is used in this work for spectral analysis and 
to compute transfer functions, i.e., to deconvolve the pulse 
generator output from the measured field waveforms. Since 
the input waveform has no frequency content above 1500 MHz, 
the transfer functions become noisy and meaningless above 
this frequency. For this reason, the transfer functions 
are low pass filtered there. 
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Prony Analysis 

The introduction of the singularity expansion method 
(SEM) [12] as a procedure for representing and calculating 
the transient electromagnetic response of antennas a^d scat- 
terers has had tremendous impact on the field of electro- 
magnetics. Its development was stimulated by observation of 
the transient responses of various complicated scatterers, 
especially aircraft. The impulse response waveforms of these 
objects are dominated by a few exponentially damped sinusoids. 

The Laplace transform of each damped sinusoid contains a 
pair of complex conjugate poles in the complex frequency plane. 
The poles represent the damped sinusoids and are intrinsic to 
the object; they are called the natural frequencies of the 
object. Under broadband excitation/ the object will have a 
large response at frequencies near its poles. The charge/ cur- 
rent/ or field distribution associated with each natural fre- 
quency is called a natural mode. The natural frequencies and 
mode distributions are characteristic of the object alone and 
are not functions of the input (incident field or current) . 

The amplitude coefficients of the natural modes of course do 
depend on the input. In the case of lightning striking an 
aircraft/ the lightning channel and the aircraft constitute 
the object. Thus a change in channel attachment points would 
change the natural frequencies. The lightning current is the 
input and determines the degree of excitation of the various 


modes . 



33 


Prony analysis [6-8] is a technique for extracting the 
poles (natural freq’'iencies) and residues (containing mode in- 
formation) of an object directly from its time-domain response. 
The method involves the use of Prony 's algorithm for fitting 
a sum of complex exponentials to a curve and some other tech- 
niques to ensure that the results are physical. Once the poles 
and residues have been obtained, the impulse response or the 
frequency-domain transfer function could be determined without 
the use of a Fourier transform. 

For a continuous signal I(t), the sum may be written 
as 


Kt) 


N 

= I 

m=l 


m 


m 


where the s_'s are the poles and the A 's the residues, both 
m m 

complex. For a discrete (or sampled) signal Kt^^), 




N 

= I 

m=l 


S nAt 
m 


m 


(3.1) 


Here, At is the sampling interval of the waveform which 
contains 2N points. N is the number of terms in the ex- 
pansion, generally referred to as the number of degrees of 
freedom in the waveform. The waveform is chosen to have 
2N points for a specific reason which will become apparent 


later. 
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For the type of transient waveforms discussed earlier 
which tend toward zero as time increases, the real parts 
of the poles will be negative. Thus, the expansion is 
seen to be a series of damped sinusoids. Also, the re- 
sulting waveform must be real, and therefore the residues 
will occur in complex conjugate pairs as will the poles. 
Prony's algorithm for computing the poles and residues for 
a given waveform is described in the following. For later 
convenience, write Eq. (3.1) as 


N S^At „ N ^ 

= I ^ )^ = j; Z" , for n = 0,1,...,2N-1 

"" m=l m=l “ 

(3.2) 

Now, consider a polynomial of the form 


(Z-Z^) (Z-Z 2 ) (2-Zj^) = 0 , 

where the Z^^'s are those of Eq. (3-2). The polynomial may 
be written as 


or 

a + a,Z + a^Z^ + ... + a.^Z^ = 0 (3.3) 

o 1 2 N 


The roots of this polynomial will be related to the poles 
by 



m 



z 


m 


(3.4) 
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Now, the a's must be found, and then the poles will be given 
by the zeros of Eq. C3.3). For the discrete case, Bq. (3.1) 
must hold for each of the 2N points in the waveform. There- 
fore 
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+ 

*2 
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+ 

> 

to 
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to 

+ 

^2 

= A^ZJ 

+ 

^2^2 

+ 




(3.5) 


I “A + A + 

^ 2 N “ 1 “ ^ 1^1 ^ ^ 2^2 ^ 




Multiplying the first equation in Eq. (3 .5) by C!^, the 

second by a^, and so on until multiplying the Nth by 

and adding the resulting N+1 equations gives the following: 



= “0*1 

+ a^Aj 

+ . . . + 


“ 1^1 

= “1*1^1 

+ Oj^AjZj 

+ . . . + 

^l\^N 

“2^2 

= 

+ a^AjZ^ 

+ ... + 

°^ 2 \^N 


N N N 

S^N " ^N^l^l ‘^N^2^2 + ..• + “AjZjj 


N 


I 

p=0 


P P 


0 


The right hand side is zero due to Eq. (3 .3) . Similarly, 
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multiplying the second equation in Eq. C3.5) by the third 

iili 

by and the (N+1) by and adding the resulting 

equations gives 


*^0^1 

= -X^Aj^Zj. 

+ 

“ 0 ^ 2 ^2 

+ 

^1^2 

= “lAj^zJ 

+ 

°‘l^2^2 

+ 

“ 2^3 

= a2A^zJ 

+ 


+ 




“N^N+f 








N 


p=0 


^p^P+1 


0 


The right hand side is again zero due to Eq. (3.3). Con- 

iili 

tinuing this procedure N times, yields on the N 


“o^N-l 

“l^N 

“2^N+1 








N-1 
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N+1 
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+ 

+ 


^^0^2^ 2 


^^1^222 


°^2^2^2 


N-1 
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N+1 
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+ 
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N-1 
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^ ^ ^ 2N-1^ . „ 2N-1, 

%^2N-1 “n^I^I '^N^2^2 


“n\^n 


2N-1 
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These N different* equations may be combined as 
N 

I ot I ,, = 0 / for k — 0 f 1 f t , , f (N~1) 
p=0 P 

If we let a.. = 1, then 
N 

N-1 

plo ”^N+k ' ^ 0,1, . . . , (N-1) 

The above equation represents N equations and N unknowns, 
the a's. These equations are written in matrix form as 



h 

^2 

T 

* • * -^N-1 

a 

o 

“^N 


^2 
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h 
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• • • ^N+1 

«2 

“^N+2 

• 
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• 

• 

• 

« 

1 — i 
1 

• !3 

H 

} 


^N+1 

^2N-2 

• 

“n-1 

• 

2N-1 

/ *5 £ 


(3.6) 


For the case considered here where the number of 
points in the waveform is equal to twice the number of 
degree’ of freedom, Eqs. (3.6) may be solved directly 
for the a‘s. Then, the roots of the N^ order polynomial 
in Eg. (3.3) yield the poles using Eq. (3.4). Using the 
poles and substituting into the first N equations of Eq, 
(3.5), the residues can be found. If the number of wave- 
form points is greater than twice the nun\ber of degrees of 
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freedom, Eg, (3. 6) will have more equations than unknowns 
and the method of least squares may be used to solve them 
approximately . 

In several ways the Prony analysis is an extremely 
useful numerical tool. If, for example, a given waveform 
is known to contain 10 poles, only 20 points of the wave- 
form are required in order to determine its poles and 
residues. The waveform can then be extrapolated forward 
in time using Eq. (3.1). Also, the data has been compressed 
so that a waveform with hundreds of sample points can be 
saved by storing only a small number of poles and residues. 

Prony analysis does suffer from some disadvantages. 

Aliasing can occur in the same manner as with the 
DPT if the sampling rate is not at least twice the fre- 
quency of the highest pole. It is also difficult to know 
how many degrees of freedom to assume for the data, some- 
times called the Prony order. If the order is assumed too 
small, the resulting poles will be incorrect. If the order 
is too large, the correct poles will be given along with 
some number of other poles which depends on how high the 
assumed order is. The presence of these extra poles will 
not greatly affect the desired poles. Our procedure for 
determining the correct order is to begin with a small order 
and increase it until consistent results are obtained when 
the order is changed. 
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The frequency resolution of the Prony analysis is 
similar to that described earlier for the FFT. The resolu- 
tion depends on the total observation time which is analyzed. 
Therefore, the maximum amount of data should be used while 
still keeping a short enough sampling interval to prevent 
aliasing [9] . 

Once the poles are calculated, they are examined for 
any which have positive real parts. These are obviously 
not physical and are removed before any residues are calcu- 
lated. Generally, these result from assuming too large an 
order and their residues are several orders of magnitude 
below the desired pole residues. In order to determine 
which of the remaining poles represent physical natural 
frequencies of the object, a technique called sliding 
window Prony [9-11] is used. The unwanted poles, called 
curve fitting poles, are put in by the algorithm just to 
fit the curve and are due primarily to noise. The sliding 
window technique involves doing several successive Prony 
analyses on a single waveform. After each, the window of 
2N points required by the analysis is moved by one sample 
point, This is normally done 3 to 10 times and the poles 
are plotted in the complex plane. The poles which represent 
physical natural frequencies will tend to cluster and the 
curve fitting poles will scatter and move around in the 
complex plane as the window is shifted. The natural fre- 
quency is taken to be the average value of the poles in 


the cluster. 
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If the input waveform from the pulse generator contains 
any poles, these will show up in the Prony analysis of the 
field measurement waveform. To avoid this, the transfer 
function is used to apply a mathematical input, containing 
known pole locations, to the system. The resulting output 
waveform is free of unknown pole contributions from the 
input. A square, pulse input is used in this work, which 
has poles only at zero and minus infinity. 

The transfer function could be obtained from the poles 
and residues by using the following expression; 


H{jw) 
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m=l 
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where the poles s_ have been written in terms of their real 

m 

and imaginary parts as 


= 
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The constant C gives rise to a delta function at t = 0. 


Prony Analysis Code 

The Prony analysis code written in Fortran is given in 
Appendix II. A flow chart of the code is shown in Figure 
3.1. Available Fortran subroutines are used to solve sets 
of linear equations and to find the zeros of a polynomial 
[24]. Due to computer limitations, a maximum order of 36 
can be handled, which corresponds to 36 poles and residues. 
After the poles and residues have been computed from each 



POOS QOal# 



Figure 3.1. Flow chart of Prony code 
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window, the waveform is reconstructed using Eq. (3.1). The 
RMS error between the original and reconstructed waveforms 
is calculated to give an indication of how well the curve has 
been fitted. Errors of up to 10% are considered acceptable, 
but 1 to 2% is typical. Detailed operator instructions for 
running the code are given elsewhere [20] . 

Stick Modeling 

In an attempt to predict the natural resonances of the 
F-106B for comparison with the Prony results, a stick-model 
technique is used [25]. The fuselage, wings, and tail of 
the aircraft are modeled using thin, perfectly conducting 
sticks. This works fairly well for normal aircraft, but the 
delta-wing does not lend itself to modeling by a single stick. 
Nevertheless, useful results may be obtained from the method. 
The model is in free space, with an electromagnetic wave polar- 
ized along the fuselage incident upon it. To first order, the 
current on each stick can be written as a directly induced 
component plus a sine and cosine component. An expression is 
written for each stick in the following form: 

"n<«> = "n -in[k(5-Vl + =°=[k(5-VJ 

The and the must be chosen to satisfy the boundary con- 
ditions at the junctions and at the ends. At the ends, the 
current must vanish. At the junctions, the current must be 
conserved and the charge density must be continuous. After 
applying these conditions, a matrix set of equations is ob- 
tained for the S and C . The natural frequencies are found 

n n ^ 
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by setting the incident field * 0 and examining these equa- 
tions for nontrivial solutions. The result is a transcendental 
equation whose zeros give the frequencies of the natural reso- 
nances. For a more detailed explanation, the reader is referred 
to reference [25] . Results from this method will be given in a 
later chapter. 



CHAPTER IV 


CIRCULAR CYLINDERS 

Prior to aprvlying the techniques described in the previous 
two Chapters to the P-106B models it was desired to test them 
to see if they give reasonable results on simpler objects. The 
natural frequencieii for isolated thin cylinders have been calcu- 
lated by Tesche [26] . If a simple circular cylinder with end 
caps were connected in place of the aircraft model, the results 
could be compared with those of Tesche for isolated cylinders. 

In this way, if natural frequencies related in a reasonable way 
to Tesche 's were found, the experimental procedure would be 
verified and at the same time valuable information obtained on 
the effect of wire attachments on natural frequencies. 

The experiments consisted of making B-dot and D-dot measure- 
ments on two different cylinders, one 2 inches and the other 6 
inches in diameter. Each cylinder is 3 feet long. The sensors 
were placed exactly in the center (lengthwise) of the cylinders. 

The reason for this placement is described as follows; The situ- 
ation is similar to that where the cylinder is isolated and a 
plane electromagnetic wave is incident upon it polarized with 
the electric field in the longitudinal direction. The modes 
with current flow in the longitudinal direction will be excited. 
These modes must have current minima or nodes and voltage maxi- 
ma or antinodes at each end. These modes are shown in Figure 
4.1. Prom this figure, it is seen that the B-dot sensor is located 
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Figure 4.1. First four cylinder modes. 
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exactly on a current node for the even modes. Similarly, the 
D-dot sensor is located on a voltage node for the odd modes. 
Consequently, the B-dot sensor will detect only the odd modes 
and the D-dot will detect only the even modes. Each wave- 
form will therefore contain only half of the natural fre- 
quencies and thus be simpler to Prony analyse. 

The measured B-dot and D-dot waveforms for the thin 
(2 in) and the thick C6 in) cylinders, with wires, are given 
in Figures 4.2 and 4.3, respectively. The input waveform 
from the pulse generator is shown in Figure 4.4. The meas- 
ured waveforms were integrated to obtain the fields and the 
input waveform deconvolved from them to produce the transfer 
function magnitudes shown in Figures 4.5 and 4.6, The cor- 
responding phases for each are shown in Figures 4.7 and 4.8. 
The peaks in the magnitude spectra indicate the resonant 
frequencies of the cylinders. The values of the first 4 
resonant frequencies for each cylinder are listed in Table 
4.1. Due to the bandwidth limitivtion of the B-dot sensor, 
only the first 4 can be reliably examined. Note the missing 
modes in the spectra for each sensor as described earlier. 

The lowest resonance at 156 MHz corresponds to the half-wave- 
length resonance of the cylinder. The cylinder is 3 feet 
(0.91 m) long, so the frequency having twice this as its 
wavelength is 

c _ 3-10^ 
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Figure 4.2. 


Measured waveforms from the thin cylinder. 
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Figure 4.3. Measured waveforms from the thick cylinder 
























TABLE 4.1. RESONANCES OF CIRCULAR CYLINDERS 


MODES PRESENT 
B 
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The frequency resolution for the spectra is 39 MHz, so 
quantization error is large enough to account for errors 
in spectral estimation. Note that for the thin cylinder, 
all higher order resonances are harmonics of the lowest 
resonance. This is not true for the thick cylinder. If 
the thin cylinder were infinitely thin, it would have harmonic- 
ally related resonances at A/2, A, 3 A/2, etc. The thick 
cylinder has more end capacitance than does the thin one 
which makes it look electrically longer and lowers the 
resonant frequencies. The Prony analysis will later show 
that the first resonant frequency of the two cylinders 
is not the same and that the thin cylinder also has its 
resonant frequencies affected by its finite thickness. 

Using the transfer functions, the square pulse was 
applied to produce output waveforms for use in the Prony 
analysis. These are shown in Figures 4.9 and 4.10 for the 
thin and thick cylinders, respectively. The waveforms from 
the thin cylinder show very sharp reflections from the ends 
as compared to those from the thick cylinder. The large 
ends cause greater dispersion of the waveform on the thick 
cylinder. The reflections die out within 15 ns, indicating 
that by then, all energy is removed from the cylinder through 
radiation and by flowing onto the wires. Note that in the 
D waveform for the thick cylinder, the effect of the end 





le-^ B (T) 



lE-9 D £C/M=) 



Figure 4.9. Field waveforms for the thin cylinder. 
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Figure 4.10. Field waveforms for the thick cylinder 
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capacitance can be readily seen as an exponential discharge 
toward zero. The cylinder and wires are like a parallel RC 
circuit with an impulse function input. The capacitor is 
charged instantaneously by the impulse and then discharges 
afterward. 

Prony analysis was next applied to the waveforms in 
Figures 4.9 and 4.10. Ten windows were used on each wave- 
form in the sliding window method. Initially, an order of 
10 was assumed and then the order was increased until con- 
sistent results were obtained. This occurred at an order 
of 20 and above. The results shown were obtained using an 
order of 20. Every 4th point was taken yielding a sampling 
interval of 200 ps and a folding frequency of 2.5 GHz, well 
above the highest frequency content of the waveforms. This 
gave an observation time of 8 ns which was enough to provide 
adequate frequency resolution. 

The results from the Prony analysis are shown in 
Figure 4.11. Four different sets of natural frequencies 
are plotted in the complex frequency plane, with the first 
4 shown for each set. The natural frequencies are also 
listed in Table 4.2. The vertical axis is the imaginary 
part of the natural frequency indicating the frequency of 
the resonance. The horizontal axis is the real part corres- 
ponding to the damping rate of the resonance. Only the 
upper left quadrant is shovm, but each natural frequency 
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Figure 4.11. Prony analysis results from the circular 
cylinders. 




TABLE 4.2. NATURAL FREQUENCIES OF CIRCULAR CYLINDEitS 


Thin 
-0.114 + 
-0.168 + 
-0.205 + 
-0.225 + 


Isolated 

Thick 

jO.869 -0.125 + jO.822 

jl.806 -0.190 + jl.772 

j2.742 -0.228 + j2.635 

j3.70Q -0.265 + j3.500 


Wire Attachment 
Thin Thick 

-0.394 + j0.980 -0.273 + jO.869 

-0„407 4- jl.960 ' -0.310 + jl.750 

-0.395 + j2.840 -0.359 + j2.643 

-0.403 + j3.840 -0.370 + j3.473 
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has its corresponding complex conjugate in the lower left 
quadrant. The natural frequencies are scaled by the factor 
L/c7t as described earlier. 

The two sets of natural frequencies closest to the 
imaginary axis are for isolated cylinders of different 
thicknesses as computed by Tesche. The leftmost sets 
are those computed here for the cylinders with wires attached. 
The information on the cylinder thickness and wire sise is 
given on Figure 4.11. Tesche did not compute the natural 
frequencies for a cylinder as thick as the 6 inch one used here, 
so those plotted are for the thickest that he did consider. 

The natural frequencies for isolated cylinders lie 
in nearly a vertical line, with damping increasing slightly 
with frequency. The natural frequencies of an infinitely 
thin cylinder would lie on the imaginary axis at 1^0, 2.0, etc. 
Therefore, the thin cylinder is closer to the imaginary 

s 

axis and more nearly vertical. As the cylinder is made 
thicker, the resonant frequencies decrease due to the capa- 
citance of the ends which make it look electrically longer. 

Tne damping also increases since the thick cylinder radi- 
ates more readily than does the thin one. The natural 
frequencies thus move down and to the left with increas- 
ing cylinder thickness. 

In the isolated case, only radiation and small ohmic 
losses remove energy from the cylinder. The addition of 
wires adds another energy loss mechanism since current can 
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flow off the cylinder onto the wires. On the thin cylinder 
with wires, the cylinder to wir, --atio is less than for the 
thick cylinder and thus the impedauje change at the junction 
is smaller. Consequently, the wires will carry off propor- 
tionately more current from the thin cylinder than from the 
thick one. This is reflected in the real parts of the natural 
frequencies, which indicate considerably more damping for the 
thin cylinder than for the thick one. The resonant frequencies 
of the thin cylinder are almost harmonically related and in- 
creased above those of the same size isolated cylinder. This 
is caused by the wires which change the geometry of the ends, 
reducing their capacitance and making the cylinder electrically 
shorter than when isolated. The same would be true for the 
thick cylinder, but is not apparent from Figure 4.11 since 
Tesche did not compute results for a cylinder as thick as 
this. If he had, they would lie slightly below and to the 
left of those shown. 

If the cylinder thickness were reduced until, in the 
limit, it became the same size as the wires, the resonances 
would vanish. This corresponds to the real part of the 
natural frequencies moving to minus infinity. The imaginary 
parts would at the same time be asymptotically approaching 
1.0, 2.0, etc. On the other hand, increasing the cylinder 
to wire ratio approaches t^e isolated case since the v/ire 
will have less and less cur^''e. - transferred to it. 
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It is interesting to compare the motion of the natural 
frequencies as the cylinder is made thicker for the isolated 
and wire-attached cases. In the isolated case, they move 
down and to the left indicating a lower resonant frequency 
and greater damping due to more radiation. However, in the 
case with the wires attached, they move down and to the right 
indicating a lower resonant frequency but less damping. There- 
fore, the effect of th® wires in allowing less energy to flow 
off for a larger cylinder to wire ratio is overcoming the greater 
radiation of the thicker cylinder. 

These results indicate that our entire method (measure- 
ment plus Prony analysis) is working well. The natural fre- 
quencies computed for the cylinders with wires differ from 
those for isolated cylinders in a logical manner. The wires 
increase the damping of each resonance by a factor of from 
1.25 to 4 since they remove energy from the cylinders. They 
also increase the frequencies of the resonances slightly since 
they make the cylinders electrically shorter. 



CHAPTER V 


F-106B SCALE MODEL 
Time Domain Reflectometry 

TDR was used to study the model in two stacres. First# 
the model was tested without its wings and tail leaving 
only the cylindrical fuselage. This was followed by a test 
with the model intact. The tests were conducted using the 
same experimental apparatus which was shown in Figure 2 . 1 , 
but with the TDR connected to the pulse generator cable. 

These tests help in understanding the effects of the wings 
and tail on the reflections which occur on the model. 

The TDR oscilloc -«ms for the fuselage alone are shown 
in Figure 5.1. The reflections from the various portions 
of the fuselage are labeled in each of the TDR traces. The 
initial rise of the trace shows a reflection coefficient of 
0.72, corresponding to an impedance of 307 ohms at the bot- 
tom of the lower wire. The gradual upward slope of the 
initial 10 feet of wire indicates a continuously increasing 
impedance (higher p) as the wave travels away from the ground 
plane. The cylinder, being a thicker conductor, represents 
a lower impedance than the wire as shown by the decrease in 
p at the cylinder bottom. The capacitance of the end cap 
is seen as an initially lower impedance then a quick exponen- 
tial charge up. 

The TDR oscillograms for the entire F-106B model are 
shown in Figure 5.2. Noted on both oscillograms are the 
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a) Overall trace 

Ground plane 
Cylinder bottom 
Cylinder top 

Cylinder bottom - 2nd pass 
Cylinder top - 2nd pass 
Overall top of apparatus 


1 2 
b) Fuselage detail 

1. Cylinder bottom 

2 . Cyl i nder top 


Figure V-1 . TDR oscillograms from fuselage only. 
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a) Overall trace 


b) Model detail 


1. Ground plane 

2. Model nose 

3. Tail 

4. Model nose - 2nd pass 

5. Tail - 2nd pass 

6. Overall top of apparatus 


1 . Model nose 

2. Start of wings 

3. End of wings 

4. Tail 


Figure V-2. TDR oscillograms from F-1063 model. 
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various portions of the model, including the added effects 
of the wings and tail. These traces are the same as those 
of the fuselage alone until the wings are encountered about 
2.4 ns after the nose. The wings represent a decreasing 
impedance, causing the reflection coefficient to decrease 
with distance along the model. Also, the main reflection 
from the end of the model comes from the end of the wings 
rather than from the end of the fuselage. This is shown 
by the shorter time (5.5 ns) the wave takes to traverse the 
model, as compared with 6 ns for the fuselage alone. This 
assumes that the wave travels at the speed of light on the 
model. The tail can also be seen as a short flat portion of 
the reflection coefficient curve about 6.5 ns after the nose. 
The overall reflection from the end of the model is smeared 
due to the largeness of the end including the wings, tail, 
and fuselage. 

The two most interesting results from the TDR experi- 
ments are first, that the end reflection from the aircraft 
comes more prominently from the end of the wings than from 
the. end of the fuselage. Second, the effect of the wings 
is to introduce a decreasing impedance as the wave propa- 
gates from noso to tail. This, in turn, makes the current 
flowing on the skin a function of position on the aircraft. 
The current is greater toward the tail than toward the nose. 
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Current Injection 

Transient field measurements on the scale model in- 
cluded the use of four different attachment point config- 
urations in order to produce a variety of possible excitation 
conditions. These configurations are shown in Figure 5.3, 
Configuration 1 is longitudinal excitation from nose to 
end of fuselage. Configuration 2 is similar except that 
the rear attachment point is on the port wing tip. Config- 
uration 3 is transverse excitation with current entering 
one wing tip and exiting from the other. Configuration 4 
is excitation in the remaining plane, with current entering 
the middle of the fuselage as in a swept stroke event and 
exiting from the top of the vertical stabili 2 er. The peak 
signal to RMS noise ratio for measurements on all configur- 
ations is greater than 30 dB, and usually is above 35 dB. 

Configurations 1 and 2 

The measured B-dot and D-dot waveforms for configur- 
ations 1 and 2 are shown in Figures 5.4 and 5.5, respectively. 
As in Chapter IV, these waveforms are integrated to obtain the 
fields and the input waveform is deconvolved from them to ob- 
tain the transfer function magnitudes shown in Figures 5.6 
and 5.7. The phases for the transfer functions of config- 
uration 1 only are shown in Figure 5.8. The frequency axes 
of the transfer functions have been scaled to the size of 
the actual aircraft. This is all that is required for simili- 
tude between model and aircraft. The model is nearly perfectly 
conducting, and thus conductivity scaling is not necessary. 
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The spectra from longitudinal B are very similar in 
configurations 1 and 2. Resonances are indicated at 7.5, 

13, 19.5, 24, and 35 MHz, the dominant being 7.5 MHz. Spectra 
from transverse B are also similar, both being dominated by 
a resonance at 19 MHz. Resonances also appear at 24.5, 35, 
and 41 MHz. The 19 MHz resonance decreases in 0 when the 
rear attachment point is placed on the wing tip. This in- 
dicates that the mode is associated primarily with the delta- 
wing. For D, the spectra are again very similar, showing 
resonances at 7.5, 19.5, 29, and 41 MHz. A poorly excited 
resonance is indicated at about 14 MHz in both spectra. 

Using the transfer functions, the square pulse was 
applied to produce the field waveforms shown in Figures 
5.9 and 5.10 for configurations 1 and 2, respectively. 

The longitudinal B waveforms are identical until the reflec- 
tions from the trailing edges of the wings and the rear 
fuselage occur. These are both contained in the first 
positive going pulse. Recall that negative longitudinal 
B indicates current flowing fore to aft and positive indi- 
cates current flowing in the opposite direction. The re- 
flection from the rear of the fuselage is larger in config- 
uration 2 than in configuration 1. This is because there is 
no longer a path for the current to leave the model from 
there, and thus more is reflected back toward the front. The 
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Figure 5.9. Field waveforms for model configuration 1 
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transverse B waveforms do not show readily identifiable 
reflections. They do show highly resonant behavior looking 
very much like the sum of a few damped sinusoids. The 
heavier damping for configuration 2 as compared to 1 is 
apparent from the waveforms. Recall that negative trans- 
verse B indicates current flow from starboard to port and 
positive indicates current flowing in the opposite direction. 
The D waveforms are again virtually identical until the 
reflections from the rear occur as in longitudinal B. These 
are both contained in the second positive going pulse. The 
reason for the almost constant negative reflection between 
the two positive pulses is the decreasing impedance along 
the model caused by the wings, as discussed in the TDR sec- 
tion of this chapter. 

Configuration 3 

The B-dot and D-dot waveforms for configuration 3 are 
shown in Figure 5.11. The transfer function magnitudes 
for the fields are shown in Figure 5.12. The spectrum from 
longitudinal B indicates resonances at 6, 19.5, and 29.5 MHz. 
Changes in the resonant frequencies after such a radical 
attachment point change are not surprising. The lowering 
of the lowest resonance is indicative of the fact that no 
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Figure 5.12. Transfer functior 
configuration 3, 
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wires now attach to the ends of the fuselage. This lowers 
the resonant frequency as was demonstrated in the cylinder 
analysis. The spectrum from transverse B shows the only 
recorded indication of a 16 MHz resonance. This dominates 
the lower half of the spectrum, which includes a well de- 
fined shoulder at 12 MHz. The upper half shows no dominant 
resonance except for the small peak oVt 41 MHz. For D, the 
spectrum is dominated by low frequencies. The peak at 6 
MHz is small and the broad peak from 13 to 20 MHz possibly 
includes 3 different resonances. This waveform is the 
only one which does not reach zero before the clear time 
is over. Therefore, this spectrum will be affected by 
leakage due to the windowing of the data. 

The field waveforms produced from the transfer func- 
tions are shown in Figure 5.13. The longitudinal B wave- 
form shows a sizable amount of current flowing longitudinally 
toward the nose of the model, even with transverse excitation. 
However, this is to be expected since the current will concen- 
trate itself at locations where the radius of curvature is 
small, such as the wing edge. The current thus flows toward 
the nose along the edge of the wing producing a longitudinal 
component. Otherwise, the field dies out very quickly after 
only about 7.5 ns. Transverse B shows current flowing 
initially from starboard to port, as expected from the ex- 
citation. As with longitudinal B, the field oscillations 
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Figure 5.13. 


Field waveforms for mode], configuration 3. 
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die out within about 7.5 ns. The D waveform has a very 
slow risetime, indicating that charging of the model is 
occurring. This charging will be discussed further in 
connection with configuration 4. 

Configuration 4 

The B-dot and D-dot waveforms for configuration 4 are 
shown in Figure 5.14. The transfer function magnitudes for 
the fields are shown in Figure 5.15. The spectrum from 
longitudinal B shows the 19 MHz resonance dominating, with 
weakly excited resonances at 14 and 35 IdHz, The 7.5 MHz 
resonance is not excited by this particular configuration 
since it produces bi-directional current flow on the fuse- 
lage. The 7.5 MHz resonance corresponds to current flowing 
in the same direction as in a half -wavelength antenna. The 
spectrum from transverse B shows resonances at 13, 30, and 
41 MHz. The D spectrum is heavily dominated by low fre- 
quencies and is uninteresting except for a small resonance 
at 19 MHz. 

The field waveforms produced from the transfer func- 
tions are shown in Figure 5.16. The longitudinal B wave- 
form indicates current initially flowing fore to aft super- 
imposed upon an exponentially decaying current. Oscilla- 
tions continue in the waveform until nearly 15 ns after 
the initial pulse. The transverse B waveform shows similar 
structure but v/ith less of an exponential component. Current 
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flows initially from port to starboard. The D waveform has 
a slow risetime and then an exponential decay to zero with 
little resonant structure on it. The explanation for the 
shapes of these waveforms is the following. The P-106B 
model can be thought of as being a capacitive blob placed 
on one conductor of a transmission line. The model is 
thus one plate of the capacitor and the ground plane is the 
other plate. Simple equivalent circuits for the apparatus 
and for the model and connecting wires are shown in Figures 
5.17 a and b, respectively. As the short pulse from the 
generator arrives at the model, it is charged quickly since 
the pulse resembles an impulse, which would charge it in~ 
stantaneously . The charging time corresponds to the rise- 
time of the electric flux density D. The charging current 
of a capacitor is proportional to the derivative of the 
voltage across it, so the magnetic flux density waveforms 
show current pulses during this charging time. In other 
words, the B waveforms are like the derivative of the D 
waveform. After the charging is completed, the capacitor 
discharges at a rate determined by the impedances in the 
equivalent circuit in Figure 5 .17b, This explains the pre- 
sence of the exponential decay in D and also in the B 
waveforms since they are proportional to the derivative 
of D. 

A summary of the resonances from the transfer function 
magnitude spectra of the four configurations is given in 
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a. Apparatus 



b. Model and wires (for t < 20 ns) 


Figure 5.17 


Approximate equivalent circuit for 
apparatus and model with connecting 
v;ires . 
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Table 5*1. Table 5.2 lists the tosonances in asesndine 
frequency order and indicates which waveforms contain which 
resonances.. The resonant frequencies listed are the average 
values from the different configurations. 

The resonant frequencies listed in Table 5*2 can be 
compared with those predicted using the stick modeling des- 
cribed in Chapter III. The model used was that shown in 
Figure 5.18, which also shows the stick lengths used. The 
resonances for this particular stick model are 6.8, 3.1, 

17.2, 24.5, 28.9, 34.7, and 41.6 MHz. Those correspond 
rather well with those in Table 5.2 except for the 11 and 
17.2 MHz values. The 17.2 may be the 19.5 MHz resonance 
in error, which is post>ible since all indications are that 
this resonance is due primarily to the flat plate delta- 
wing, which is not well modeled by the sticks. It could 
also bo the 16 MHz resonance slightly in error. This 
would mean that the 19.5 MHz resonance is not predicted at 
all. The 11 MHz is likely the 12.8 MHz resonance in error, 
again due to the wings. The general nature of these reso- 
nances agrees with those computed for the B-1 aircraft [ 25] . 
Specifically, the second resonance is approximately 50% above 
the first. Overall, the stick modal results tend to agree with 
the locations of the resonances as determined from the trans- 
fer functions. 

Frony analysis was next applied to all 12 field wave- 
forms from the four configurations. The order was initially 



TABLE 5.1. 


EESONANCES OF F-106B MODEL 


Configuration 1 Resonances, MHz 

Conf igur at ion 

2 Resonances 

, MHz 

®L 


°F 

^L 


°F 

7.5 

19 

7.5 

7.5 

19 

7.5 

13 

24.5 

19.5 

13 

24.5 

19.5 

19.5 

35 

29 

19.5 

35 

29 

24 

41 

41 

24 

41 

41 

35 



35 



Configuration 3 

Resonances, MHz 

Configuration 4 Resonances, MHz 

®L 

®T 

^F 

®L 


^F 

6 

12 

6 

19 

13 

19 

19.5 

16* 



30 


29.5 

41 



41 



only occurrence 



TABLE 5.2. RESONANCES OF F-106B MODEL 
Config. 1 Config. 2 Config. 3 Config. 4 


Resonance 

®L 





°p 

°L 


°F 

^L 
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X 
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X 

X 


X 


13 

X 



X 




X 


X 

16 








X 



19 

X 

X 

X 

X 

X 

X 

X 



X 

24 

X 

X 


X 

X 
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X 
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Figure 5.18. 


Stick model configuration. 



assumed to be 20 and increased until the maximum order of 
36 was reached. Consistent results were just beginning to 
be seen at this high order and so it was used for all model 
waveforms. Only three windows were used in the sliding 
window method to save computer run time, which becomes 
large when such a high order is used. It is still easy to 
distinguish the physical from the curve-f itting poles using 
only 3 windows. Every 6th point was taken in the waveforms 
yielding a sampling interval of 300 ps and a folding fre~ 
quency of 1.7 GHz, still above the highest frequency con- 
tent of the waveforms. 

The dominant natural frequencies for configurations 1 
and 2 are shown in Figure 5.19. Those for configuration 3 
and 4 are shown in Figures 5.20 and 5.21, respectively. 

The values of the natural frequencies are listed in Table 
5.3 for each configuration. These values are averages de- 
rived from the sliding window method. Each natural fre- 
quency for a given configuration was taken from v;hichever 
waveform had that particular resonance most dominant. 

Figure 5.19 shows 4 dominant natural frequencies. Note 
that as the wire is moved from the fuselage to the wing, the 
first two show a decreased damping rate. This indicates that 
the modes had current flowing off of the rear wire, which was 
stopped when the wire was removed from the fuselage. There- 
fore, these modes both involve current flow on the rear of 
the fuselage. The third shows the opposite effect with 









TABLE 5.3. NATURAL FREQUENCIES OF F-106B MODEL 


Configuration 1 


Natural frequency 

f ,MHz 

-.25 + jO.93 

8.1 

-.34 + jl.50 

13.1 

-.29 + j2.09 

18.2 

-.28 + j2.80 

24.4 

Configuration 3 
Natural frequency 

f ,MHz 

-.22 + j0.80 

7.0 

-.26 + jl.47 

12.8 

-.21 + jl.84 

16.1* 

-.36 + j2.23 

19.5 

-.42 -f j3.35 

29.3 


Configiiration 2 


Natural frequency f,MHz 

-.23 + jO.88 7.7 

-.28 + jl.45 12.6 

-.35 + j2.09 18.2 

-.27 + j2.98 25.9 

Configuration 4 
Natural frequency f,MHz 

-.22 + jl.55 13.5 

-.33 + j2.15 18.8 

-.49 + j4.00 34.9 


* only occurrence 
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heavier damping. This mode is due primarily to the delta- 
wing and so the current flows off when the wire is attached 
to the wing tip. The fourth natural frequency was poorly 
excited compared to the first three and shows a change in 
frequency rather than damping. 

Figure 5.20 for configuration 3 includes 5 dominant 
natural frequencies. The one which appears out of line 
with the others is the one at 16 MHz in its only appearance. 
The other configurations do not excite this mode. 

Figure 5.21 for configuration 4 includes only 3 domi- 
nate natural frequencies. Note again that the lowest one 
is not present indicating that it was not excited in this 
particular configuration. 

The values of the natural frequencies in T3i,.le 5.3 
correspond well to the transfer functions. The damping 
rates compare favorably with the width of the peaks and 
the resonant frequencies compare well with the centers 
of the peaks except in configuration 3, where the Prony 
analysis does not show the low 6 MHz resonance but puts 
it at 7 MHz. The frequency has indeed been lowered in 
this configuration, but the transfer function spectrum 
has too much quantization error to obtain the resonant 
frequency accurately. The Prony analysis is the more 
accurate spectrfVl estimator of the two methods. The 41 
MHz resonance was also extracted by the Prony analysis, 
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but it was not Included because it is above the bandwidth o£ 
the B~dot sensor and is probably unreliable. 

The Prony analysis has yielded natural frequencies which 
agree well with the resonances shown in the transfer function 
spectra. Comparisons of these results with computations on 
direct and nearby strike data in the next chapter will shed 
new light on the nature of the lightning channel. 



CHAPTER VI 


COMPARISONS WITH IN-FLIGHT DATA 


Overview 

A discussion of the in-flight lightning strike data 
obtained in the 19Sn and 1981 campaigns is given in refer- 
ence [14] . A total of 20 strikes to the aircraft were re- 
ported along with several nearby strikes. Observations dis- 
cussed here are for magnetic and electric fields measured 
with the longitudinal B-dot sensor over the st$irboard wing 
and the nose D-dot sensor. The bandwidth of these measure- 
ments is 50 MHz with a sampling interval of 10 ns. Lightn- 
ing pulse lengths ranging from 25 ns to 7 us have been ob- 
served. 

The strikes seen thus far have had relatively small 
peak currents, typically hundreds of amperes. On one 
occasion, however, an I sensor on the noseboom recorded a 
large stroke of 14 kA peak current shortly after the B-dot 
sensor had triggered on a very small pulse. Thus, it appears 
that streamering on the aircraft may be occurring which trig- 
gers the derivative sensors and causes subsequent large-current 
strokes to be missed. Even so, these short, small exci- 
tations contain high frequencies v;hich excite aircraft reso- 
nances . 
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Direct Strikes 

Five examples of B-dot records from direct lightning 
strikes to the aircraft are shown in Figures 6.1 through 
6.5, along with their frequency spectra, magnitude only. 

Two D-dot records are shown in Figures 6.6 and 6.7. All 
records are from 1980 flight 38 except for the D-dot in 
Figure 6.7 from flight 18. 

In some cases, particularly in Figures 6.1 and 6.2, 
the spectra indicate that aliasing is present since the 
magnitudes do not approach zero at the 50 MHz folding fre- 
quency but are actually increasing. Although the sensor 
outputs are low-pass filtered at 50 MHz, the filter attenua- 
tion is rather gradual and some higher frequencies are allowed 
to pass through to the digitizers. Probably this affects 
only th<?' higher frequency range of the spectra between 25 
and 50 MHz. 

Comparisons of the spectral peaks with those for the 
model show excellent agreement except for the 19 MHz reso- 
nance, which appears consistently at 21 MHz in the in-flight 
data. This resonance is due to the delta-wing and the discrep- 
ancy might be explained by the fact that the elevens on the rear 
of the wing are somewhat electrically isolated from the rest 
of the aircraft. They make up about 10% of the wing and aft 
fuselage surface area which accounts for the frequency 
increase on the F-106B. Two examples of a 16 to 18 MHz 
resonance are seen. One is the small 16 MHz peak in the 
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Figure 6.1. B-dot record from direct strike 80-038-01, 
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Figure 6.5. B-dot record from direct strike 80-038-04, 
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D-dot record from direct strike 80-038-04. 










109 



r-'i 

U 



f 

n 


B-dot spectrum of Figure 6.5 and the other is the 17*5 MHz 
peak in the D-dot spectrum of Figure 6.7. The higher fre- 
quency peaks in the various spectra also show the higher 
resonances seen on the model. A summary of the resonant 
frequencies indicated in the direct-strike data is given 
in Table 6.1. All model resonances are seen except for 
the 24 MHz resonance where only a small shoulder is indi- 
cated in 3 spectra. B-dot 38-4 and D-dot 38-4 have spectral 
peaks below and above the 29 MHz resonance, respectively. 

Prony analysis has been applied to all the above direct- 
strike records. The methods used were identical to those 
used in previous chapters. A difficulty to be expected in 
this is that the lightning, which is the input to the air- 
craft-channel system, will contribute some natural frequencies 
to the measured waveforms. These must be fitted by the Prony 
analysis along with those from the system, presenting added 
difficulty to the extraction, of accurate natural frequencies 
of the system. 

The results from these Prony analyses are given in 
Table 6.2. The natural frequencies from flight 38 B-dot 
records 1 and 5 are almost identical. This was expected due 
to the similarity of the spectra from the two. Of the three 
natural frequencies from record 3A, the second agrees well with 
the corresponding one in records 1 and 5. The first has greater 


TABLE 6.1. F-106B RESONANCES FROM IN-FLIGHT DIRECT LIGHTNING STRIKES 


Model B-dot D-dot 


Resonance 

38-1 

38-5 

38-3A 

38-3B 

38-4 

38-4 

18-1 

7 

6.2 

6.2 


7,0 

7.0 


8.0 

13 

13.0 

13,0 

12.0 

12.5 




16 





16.0 


17.5 

19 

21.0 

21.0 

22,2 

22.0 

21.0 

22.0 


24 

* 

* 


* 




29 





27.0 

31.0 

29.0 

35 



34.0 




35,0 

41 




40.5 


40.5 



* shoulder 
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TABLE 6.2. NATURAL FREQUENCIES FROM DIRECT-STRIKE IN-FLIGHT DATA 


B-dot 


38-1 



38-5 



Natural frequency 

f ,MHz . 


Natural frequency 

f ,MHz 


-.215 + jl.50 

13.1 


-.21 f jl.54 

13.4 

o P. 

-.295 + j2.50 

21.8 


-.31 + j2.48 

21.7 

tJ 

o 






c , = 

38-3A 



38-3B 



Natural frequency 

f ,MHz 


Natural frequency 

f ,MHz 

^ ' 

-.26 + jl.34 

11.7 


-.29 + jO.73 

6.3 

^ ■ 
»C' ■ 

-.32 + j2,48 

21.6 


-.26 + jl.28 

11.1 


-.22 + j4.00 

34.9 

38-4 

-.37 + j2.53 

22.1 



Natural frequency 

f ,MHz 




-.23 

+ j0.90 

7.8 




-.27 

+ jl. 61 

14.1 




-.30 

+ 

LJ. 

• 

O 

to 

35.1 




-.47 

+ j4.61 

40.2 





.D-dot 




38-4 



18-1 



Natural frequency 

f ,MHz 


Natural frequency 

f ,MHz 


-.29 + j2.54 

22.1 


-.08 + jO.86 

7.7 



-.13 + j2.03 17.7 


m 
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damping and lov/er frequency and the third was not present 
in either previous record. Record 3B is the first to show 
the aircraft's lowest natural frequency. It has greater 
damping than that seen before for the lowest resonance in 
either in-flight or model data. The third is also highly 
damped compared to other in-flight results. Only the first 
natural frequency from record 4 corresponds well to a peak 
in its spectrum. This spectrum is unusual since it is 
almost smooth except for a dip at 11 MHz, The other natural 
frequencies at 14, 35, and 40 MHz do not correspond well 
to peaks in the spectrum but do correspond fairly well with 
aircraft resonances, however. The domibant natural frequency 
in D-dot record 4 agrees with previous results. 

D-dot record 1 from flight 18 shows the lowest damping 
recorded in the in-flight data. It is also unusual in that 
it contains resonant oscillations with a step occurring after 
they begib. The reason for the low damping is unclear unless 
it is indicative of more than one lightning pulse. 

Nearby Strikes 

Results from nearby strikes can be used to determine the 
natural frequencies of the F-106B for the case of no channel 
attachment. This is the limiting case where the damping is 
the lightest. Again, the lightning itself is the system 
input for these records and problems may occur in the accur- 
ate extraction of the aircraft resonances . 
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B-dot and D-dot records from a nearby lightning strike 
are shown in Figures 6.8 and 6.9, respectively, along with 
their corresponding magnitude spectra. The B-dot record 

ii 

utilized only 8 levels of quantization in the largest peak 
whereas the D-dot record used 31. The B-dot will thus suf- 
fer more from quantization error. Aliasing is also present 
in both magnitude spectra. The first four aircraft reso- 
nances are shown at 7.3, 11.9, 15,4, and 21 MHz, the 15.4 
MHz being only a shoulder in each spectrum. These are in 
agreement with those obtained from the model and previous 
in-flight data. D-dot waveforms and spectra from two other 
nearby strikes are shown in Figures 6.10 and 6.11. These 
waveforms are unusual because of their unipolar shape, 
which indicates that D increased to some value and remained 
there. The spectra show the first aircraft resonance at 
about 7 MHz. All the nearby strike records are 1981. 

The first two records are from the same lig*.. . event, 
so the natural frequencies should ag:?fffie. They do agree ex- 
cept for the damping on the 15.3 MHz resonance and for the 
frequency on the 21 MHz resonance. Neither of these reso- 
nances are as well excited in the D-dot record as they are 

Ij in B-dot, so the B-dot results are the most reliable. This 

is supported by the fact that the 21 MHz resonance in B-dot 

ii 

11. agrees with all previous in-f light results. 
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Figure 6,8. B-dot record fro.i« nearby strike !/l-026-lQ, 




FPEaUB>CY (hHZJ 


Figure 6.9. D-dot record from nearby strike 81-026-10, 
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Figure 6.11. D-dot record from nearby strike 
81-026-07. 
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Of greater interest are the results from Prony analysis 
since the damping for each resonance should be at its minimum 
value. These results are listed in Table 6.3. The dampings 
of the resonances are indeed lighter than those obtained pre- 
viously from either the model or direct-strike in-flight data, 
indicating that the dampings are representative of the isolated 
aircraft. The only exception is the unusual D-dot record 
from flight 18 . This is an anomalous case which cannot be 
explained using these techniques. In regard to the natural 
frequencies extracted from records 81-026-03 and 07 it should 
be pointed out that, because of the coarse quantization in- 
herent in these waveforms, the Prony anaJ.ysis is somewhat 
approximate. The RMS error between the waveforms and their 
Prony reconstructions ranged up to 20%, so there is some un- 
certainty in the values of the natural frequencies. 

The heaviest damping of the resonances should occur on, 
the model, since it has fairly large, highly conducting wires. 
This does occur except for the first and third natural fre- 
quencies from B-dot record 38-3B. Perhaps multiple channels 
caused the high damping in this case. 

The natural frequencies from flight 38 direct-strikes 
and from nearby-strikes are plotted in the scaled complex 
plane in Figure 6.12. Also shown are the results from model 
configuration 2 since it is most similar to the actual attach- 
ment configuration for these strikes. The third natural fre- 
quency from the nearby-strike records is not shown since 
there is no corresponding frequency in the direct- strike 
results. The natural frequencies from the lone D-dot record 


TABLE 6.3. 


F-106B NATURAL FREQUENCIES FROM NEARBY LIGHTNING STRIKES 


B-dot 26—10 
Natural frequency 

f ,MHz 

D-dot 26-10 
Natural frequency 

f,bHIz 

-.19 + jO.84 

7.4 

-.14 + jO.84 

7.3 

-.16 + jl.35. 

11.8 

-.15 + jl.36 

11.9 

-.13 + jl.76 

15.4 

-.24 + jl.76 

15.3 

-.21 + j2.41 

21.1 

-.17 + j2.19 

19.1 

D-dot 26—3 
Natural frequency 

f ,MHz 

D-dot 26-7 
Natural frequency 

f ,MHz 

- . 14 + j 0 . 79 

6.9 

-.16 + j0.71 

6.2 
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H* 

\o 


ORIGIi'Wt. PAGE (3 
OF POOR QUALITY 



0{; POOR QUALITY 


In-flight Direct 
3-dot 

□ 80-038-01, 80-038-05 

A 80-038-03A 
0 80-038-03B 

+ 80-038-04 

D-dot 

V 80-038-04 


• In-flight Nearby 
(average) 

X Model Configura- 
ation 2 






121 


from flight 18 are plotted in Figure 6.13 along with model 
configuration 3 since it was the only one which showed a 
16 MHz resonance. 

The general tendency is for the in-flight results to 
show lighter damping than the model. This indicates that 
possibly the streamer theory discussed earlier is true, 
meaning that the aircraft is almost isolated. Another 
possibility is that a second attachment pcrlnt, or a detach- 
ment point, has not occurred since the leader probably 
travels only a few cm in the time of most events, about 1 
us. Therefore, only one attachment is present and the 


4. csaOuciiiwcso cixc: uiwJ.C3 Xx^iluJuj^ dauip^\a • r% uUXa-vi 

possibility is that the channel, being lossy, is inducing 
higher Q in the resonances since the current cannot flow off 
the aircraft as readily as on the highly conducting wires 


on the model. 
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CHAPTER VII 
CONCLUSIONS 

The laboratory modeling and analysis techniques des- 
cribed in this report constitute a high signal-to-noise ratio 
method for studying the external response of an aircraft to 
injected current pulses. This has been shown to be a useful 
approach to modeling direct lightning strikes to the aircraft. 
The idea is not to try to duplicate in the laboratory the 
waveforms observed in a direct strike, but rather to deter- 
mine the resonances or natural frequencies of both the labora- 
tory and in-flight waveforms and then to use the comparison 
between them to help in the interpretation of the in-flight 
results . 

As an initial test of the method, application to circular 
cylinders has provided new information on the effect of wire 
attachments on scatterers as compared to their isolated be* 
havior. The wires cause the damping of each resonance to 
increase since they remove energy from the cylinder, adding 
to that removed by radiation. They also increase the fre- 
quencies of the resonances slightly since they make the 
cylinder electrically shorter. 

Using the P-106B aircraft model, good agreement has been 
seen in comparisons with the rather limited amount of in-flight 
data that is available. Resonances excited by direct and 
nearby lightning strikes are substantially the same as those 
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on the model. However, detailed comparisons of model 
natural frequencies with those computed for direct-strike 
data almost always show higher Q in the in-flight case. This 
is indicative of either a lightning channel with a higher 
impedance than the wires on the model, only one channel attach- 
ment point, or short streamers instead of a long channel. 

More in-flight data may help to sort this out. For the nearby 
strike data, where the aircraft is isolated, the Q is the 
highest (except for one anomalous case) . The average damping 
for the lowest resonance is -0.16 as compared to -0.24 for 
the model with its highly conductinc; wires. The value of 
-0.16 for the isolated aircraft is comparable to that of 
an isolated cylinder with D/L - 0.17 but is semewhat higher 
than typical stick model results for various aircraft [25] . 
Quantization errors limit the accuracy of the Prony analysis 
on the in-flight data. 

Twelve waveforms were recorded on the model under short- 
pulse excitation conditions for the determination of the 
natural frequencies. (See figures 5.9, 5.10, 5.13, and 5.16) 
Four different wire attachment configurations v/ere used. 

The waveforms contain information on characteristics to be 
expected on the airplane such as the ratios of E to H and of 
one component of H to the other. The natural frequencies ex- 
tracted from these waveforms constitute a set of basic 
parameters for comparison with future in-flight data. (See 
Table 5.3.) 
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One area which might be addressed in future research 
is that of achieving a closer approximation to the 
lightning channel by modifying the wires on the model. This 
would involve using wires of smaller diameter and/or higher 
resistance and observing whether the natural frequencies of the 
model could be brought into agreement with those in-flight. 

If the present modeling technique were apolied to other 
types of aircraft, one would expect generally similar results. 
The similarity would probably be strongest for other aircraft 
which have the delta-wing shape of the P-106B. These include 
the shuttle, the P-16XL, various foreign-made aircraft, and 
several variable-wing aircraft such as the B-1 which are 
similar to a delta-wing when their wings are in the swept 
position. One would expect in the case of an aircraft contain- 
ing large amounts of composite materials that the resonances 
would be more highly damped due to the lov;er electrical 
conductivity. 
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APPENDIX I 


LISTING OF DATA ACQUISITION CODE 



DATA ACQUISITION CODE «*»*******^Ht*»» 


DJfIFNSJON IDATAA(5j2) , IDATAB(512) » IDATAC(5J2) r IDATAD(512) r 
IDATAE(512) r IDATAF(512) , IDATAQ(512> rIDATAH{512) , 

JDATAJ (5J2) »IDATA.H512) ,IDATAK(5J23 rRDATA(5J2) , A(3) 

BYTE NAME(IB) 

RLSB=^ .a028E-3 
CAI^L DATE(A) 

GO TO 60 
TYPE 200 
ACCEPT 510, L 
IF(L.EQ.O)GO TO 300 
TYPE 300 
ACCEPT 310, NAME 
TYPE 400 

ACCEPT 410,NSAMPL 
TYPE 500 

ACCEPT 510,IFIRST 
TYPE 710 
ACCEPT 410,IDIM 
IF( IDI0.EQ.200)SGAIN=1 .02 
IF( IDIU.EQ.i00)SGAIN=2.0Q 
IF(IDIU.EQ.50)SGAIN=4.1 
IF(IDIM,EQ.20)SGAJN=10.2 
IF( IDIV.EQ.10)SGAIN=20.25 
IF( IDI0.EQ.5)SGAIN=41.0 

IF( IFIRST.EQ.0.0R.IFIRST.EQ.1)GAIN=10.34«SGAIN 
IF( JFJR5T.EQ.2.0R. IFIRST.EQ.3)GAIN=J0.08*S6AIN 
IF(IFIRST.EQ.4)GAIN=2.01*SGAIN 
IF(IFIRST.EQ.5)GAIN=2.01»SQAIN 
IF n FIRST. EQ.6)GAIN=108.8*SGAIN 
IF( IFIRST.EQ. 7) GAIN= 105. 5«SGAIN 
TYPE 600 

ACCEPT 5J0,NCHANN 

- tj 
o 
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n n 


TYPE 700 

ACCEPT 410 ,IM 0 DE 
NCHANN=1 
IMDDE =0 
COR^RLSB/GAIN 
INC -^1 

DO 5 1 = 1,512 
IDATAA( I )=0 
IDATABl I )=0 
IDATAC(I )=0 
IDATAD(I )=0 
IDATAE(I )=0 
IDATAF(I )=0 
IDATAG(I )=0 
IDATAHC I )=0 
IDATAKI )=0 
IDATAJ<I )=0 
IDATAK( I )=0 
5 RDATA(I)= 0 . 

TYPE 850 

ACCEPT BSOrSWTIM 
TYPE 870 

ACCEPT BBOrSWLEN 
RATE= ( 40 . *SWLEN ) /SWTIM 
RC 0 NT= 1 .EB/RATE 
ICOUNT=INT(RCONT) 

RCOUNT=FLOAT{ ICOUNT) 

1 F( (RCONT-ICOUNT) .GE. 0 . 5 )RCQUNT=RC 0 UNT+ 1 . 

RATEAC= 1 . EB/RCOUNT 

SINT= 2000 . / ( (RATEAC^SWTISI) /SWLEN) 

TYPE 890 , RCOUNT , RATE , RATEAC , SI NT 

IFLG =0 

IFLOJ =0 

CALL RTS(IDATAA, 512 , rNSAMPL, IFIRST,iNCHANNr , IMODE, IFLG, IDUM) 
CAM. SETRM , J 3 , RCOUNT , IFLOl 1 


o o 

■P ^ 
O - ? 
© 

^ . 

C r.3 




•< £3 


LiJ 

H* 



CALL LWAITCIFLGrO) 

CALL SETR(-1,,,) 

IFLG^O 

CALL 8TS ( I D ATAB r 5 1 2 , , NSAHPL r I F I RST r NCHAKM . , I MODE , I FLG , I DUM > 
CAU . SETR ( 3 r J 3 , RCOUNT , I FLG J ) 

CALL LWAIT(IFLG,0) 

CALL SETR(-1,,,) 

IFLG=0 

CALL RTS(IDATAC#512r »NSAMPL» IFIRST »NCHANNr » IMODE^ IFLGrIDUM) 
CAU SETRC1,13,RCQUNT,IFLQ1) 

CALL L«AITaFLQ,0) 

CALL SETR(-1,,,) 

IFLG=0 

CALL RTSCIDATAD,512, »NSAMPL, IFIRSTrNCHANN, , IMODEr IFLQ, IDUM) 
CAU HETR (1,13, RCOUNT , I FLGl ) 

CALL LWAIT(IFLG,0) 

CALL SETR(-1 , , , ) 

IFLQ=0 

CALL RTS(IDATAE,512, rNSAMPL, IFIRST, NCHANN, , IMODE, IFLG, IDUM) 
CALL SETR( 1 , 13,RC00HT,JFLGn 
CALL LWAIT{IFLG,0) 

CALL SETR(~1,,,) 

IFLG-^0 

CALL RTS(IDATAF,512,,NSAMPL,IFIRST,NCHANN, , IMODE, IFLG, IDUM) 
CAU SETRU , 13, RCOUNT, IFLGl ) 

CALL LWAIT(IFLG,0) 

CALL SETR(-1,,,) 

IFLG-0 

CALL RTS( IDATAG,512, ,NSAMPL, IFIRST, NCHANN, , IMODE, IFLG, IDUM) 
CAU SETR ( 1 , 1 3 , RCOUNT , I Fl,G 1 ) 

CALL LWAIT(IFLG,0) 

CALL SETR(-1,,,) 

IFLG=0 

CALLRTCTIDATAH,512,,NSAMPL,IFIRST, NCHANN , , I MODE , I FLG , I DUM ) 
CAM SFTR ( 1 , 1 3 , RCOUNT , J FLG 1 ) 



CALL LWAIT(IFLGrO) 

CALL SETRC-l r , , ) 

IFLG=0 

CALL RTS( IDATAI ,512r ,NSAMPLr IFIRSTrNCHANN, r IMODE» IFLQr IDUM) 

CALL SETR< J , J3rRC0UNT, JFLGJ ) 

CALL LWAIT(IFLGiO) 

CALL SETR(-1 ,, ,) 

IFLG=0 

CALL RTS( IDATAJrSlZr »NSAMPt., IFIRSTrNCHANN, , IMODE, IFLG, IDUM) 

CALL SETR () , ) 3 , RCOIJNT , I FLG J ) 

CALL LWAIT( IFLGfO) 

CALL SETR(-l,r,) 

IFLG=0 

CALL RTS( IDATAK.SIZ, ,NSAMPLr IFIRSTrNCHANN» » IMODE, IFLG, IDUM) 

CAI L 5FTRU , J3»RC0UNT» JFLGl ) 

CALL LWAIK IFLGrO) 

CALL SETR(-1 r , r ) 

DO 10 I=1,NSAMPL 

RD ATA ( I ) = I D AT AB ( I ) + 1 DATAC ( I ) + 1 DATAD ( 11 ) + 1 DATAE ( I ) + 1 D ATAF ( I ) 

J 0 RDATA( J )=RDATA( I ) + IDATAG( I ) + IDATAH( I H+IDATAI ( I ) + IDATAJ( J ) 

X +IDATAKU) 

DO 12 I=1.NSAMPL 

12 RDATA(I)=( (RDATA(I)/10. )-2048. )*COR 

TYPF 320,NA.1F 
TYPE 330, A 
DO 15 I=lrNSAMPL,5 

J 5 TYPE 20 , RDATA ( I ) r RD ATA ( I + 1 ) , RDATA ( I +2 ) , RDATA ( I +3 ) , RDATA ( 1+4 ) 

OPEN (UNIT=2,NAME=NAMErTYPE= 'UNKNOWN '.-REC0RDSI2E=J , 

X ACCESS= 'DIRECT ',ASSOCIATEUARIABLE=INC) 

DO 37 J=J ,532 

37 NR]TE<2'INC) RDATA(I) 

CLOSE (UNIT=2, DISPOSES 'SAVE ' ) 

GO TO 50 

20 F0RMAT(1X,5E12.4) 

200 FORMAT! IX, 'SAMPLE ANOTHER WAVEFORM? '/ 3 X, 'TYPF 3 IF YrO IF NT ',$) 


50 

T3 O 

51 

S3 fs: 

o 

C Pi 

■< (f) 


»-• 

to 


300 FORMAT (IX, "ENTER FILE NAMErAS "DYO »KXXXXX . XXX" ' » $ ) 

310 FORMAT! 16A1) 

320 FORMAT ( /2GX , 1 BA 1) 

330 FORMAT { 27X, 3A4 ) 

400 FORMAI ( IX, "TAKE HON MANY SAMPLES? ",$) 

410 FORMAT! 15) 

500 FORMAT ! IX, "FIRST CHANNEL TO BE SAMPLED?", $) 

510 FORMAT! ID 

600 FORMAT! IX "HOW MANY CHANNELS TO BE SAMPLED?", $) 

700 FORMAT ! IX, "MODE?" ,$) 

710 FORMAT! IX, 'OSCILLOSCOPE MILLJMOLTS/DIM?" ,$) 

H50 FORMAT! JX, "SWEEP TJME?",$) 

870 FORM AT! IX, 'SWEEP LENGTH? ",$) 

BGO FORMAT! FI 0.2) 

BSO F0RMATI/5X, "RCOUNT= ",FIO. 2, /5X, "DESIRED SAMPLING RATE 

X F10.2," Hz ",5X, 'ACTUAL SAMPLING RATE= " ,FJ 0.2, " Hz"/ 

X 5X, 'SAMPLING INTERVAL^ ' ,FJ 0. 5, ' f»S ' ) 

CALL EXIT 
END 


a 00 
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C 

C SINGULARITY PKPANSinN 

C 

C USING 

C 

C PRONY'S METHOD 

C 

C MAX I HUM OF 36 DEGREES OF FREEDOM 

C 

C^^*-it»^-!Hf»»*ii-»^*-lt*********1Hf*^*»»**if»*******1t**it** 
COMHON/l^iTA/RDATA (512) 

COMMON/READ/TDATA ( 512 ) 

COMMON/NDATA / A VD ATA ( 72 ) 

COMMON/PARA/N,NN,NS,NSHJF 
COMMON/DAY/NAME r A 
COMMON/PRO/DELTA , DELT 

DIMENSION WW(37) ,R00TRE(3B) rR00TIM(3B) r ALMAT(3B,3B) ,ALPHA(37) 
BYTE NAME(16) , AC IB) 

CALL DATECA) 

TYPE 125 

ACCEPT 310,NRESUf1 
IF(NRESUM.EQ.O)GO TO 1 
TYPE 128 

ACCEPT 310,NEWSFT 
NEWSFT=NEWSFT+1 
GO TO BO 
1 TYPE 200 

ACCEPT 21 Or NAME ' 

TYPE 320 
ACCEPT 330 r DELT 
TYPE 325 

ACCEPT 310rNREC0N 
TYPE 205 

ACCEPT 2J0rPNAME 


S o 

Q ^ 

iO . 

•< m 


UJ 

<n 





c 

READ DATA FROM DISK 
C 

CALI-' READ (NAME) 

C 

CiHHH* INITIALIZATION **«**»»**«*»♦«««•»♦»♦*♦ 
C 

2 NEWSFT=0 
NRESUM=0 
TYPE 220 
ACCEPT 310,N 
M=N+1 
NN = 2^jN 
TYPE 300 
ACCEPT 310, NS 
TYPE 340 

ACCEPT 350,DCLEM’ 

NSAMPL=NS*NN 

IF(NSAI1PL.GT.432)Q0 TO BO 
DELTA=NS»DELT 
C 

DO SHIFTING 

c 

TYPE 100 

ACCEPT 310,NSHIF 
BO TYPE 150 

ACCEPT 310,NSHIFS 
TYPE 4B, DELTA 

4B FORMATUK. 'DELTA= ',E15.5) 

(32 NSHIFT=NHHIF-3+NEWSFT 
DO 70 J=1,NSHIFS 
CALL SETERR{10,1U 
NSHIFT=NSHIFT+1 
DO 10 1=1,452 

10 RDATAl I )=TDATA< I+NSHIFT)+DCLEO 


LO 

-J 


0Kiis:w^\L pmE B 

OF POOR QUALITY 


c 

FORM REDUCED WAMEFORM **»•»»«*# 

ADD DC LFVEL JF REOUJRED «***##* 

C 

TVPF ^5',NShiIFT 

45 FGRMAT(/1X» 'NSHIFT='rI5) 

DO 30 I^l.NN 
K = I-1 

30 A0DATA(I)=RDATA(1+(K#NS) ) 

C TYPE 47 

47 FORMAT (2BX, "AVDATA ' > 

C DO 43 1 = J rNN,5 O C 

C43 TYPE 50,A0DATA(I) ,AUDATA<I+J ) »AVDATA(I+2) ,AVDATA(I+3) 

C X,AVDATA( 1+4) 

50 I- ORMAT (1 K r 5E 3 2 . 4 ) 

C 

CALL PRONY SUBROUTINE **+*#»*»##*Jf 
C 

70 CALL PRONY(ALMAT, ALPHA, R00TRE,R00TII1,WMrN,M) 

GO TO 90 
80 TYPE 700 

90 TYPE 500 

ACCEPT 330rI00ER 
IF( lOOER.EQ.l )G0 TO 2 
TYPE BOO 

ACCEPT 310, I OVER 
IF<I0VER.EQ.1)G0 TO 1 
GO TO 900 

100 FORMATUX, 'NUMBER OF POINTS TO SHIFT(+=SHIFT LEFT?-=SHJFT 
X RIGHT)? 'rSi) 

125 FORMATC IX, 'RESUME ON PREVIOUS WAVEFORM? ',$) 

12G FORMAT! IX, 'SHIFTS TO SKIP? ',*) 

1 50 FORMAT OX.' NUMBER OF SHIFTS? ' , $ ) 

200 FORMATOX, 'ENTER DATA FILE NAME, AS ‘'DYO:XXXXXX.XXX" ',$) 

205 FORMATOX, 'ENTER POLE FILE NAME, AS "DYO:XXXXXX,XXX" ',$) 



H 

LJ 

03 


210 FORMATdBAl) 

220 FDRflAK IK r 'ORDER? ',$) 

300 FORMATUX, 'TAKE EVERY ?t-h POINT? ',$) 

:310 R]RI1AT( iri) 

320 FORMAT nXr 'WAVEFORM SAMPI.JNG INTERVAL? 'r$) 

323 FORMAT( IXr 'RECONSTRUCT AND COMPUTE RMS ERROR? ',$) 

330 FORMAT ( F20 . 3 ) 

0 FORMAT ( 1 X , ' DC LEVEL ADDED? ' 

350 FORMAT (FJ 0.5) 

500 FORMAT < IX r 'ANOTHER EXPANSION ON SAME WAVEFORM? ',$) 

GOO FORMAT! IX, 'ANOTHER EXPANSION ON DIFFERENT WAVEFORM? ',$) 

700 FORMAT! IX r 'REDUCED WAVEFORM EXTENDS TOO LONG IN TIME') 

.300 CALL EXIT 
END 

SUBROUTINE READ(NAME) 

ROUTINE TO READ DATA FROM DISK ♦***»«»*#»***iHH!*****#* 
CQMM0N/READ/TDATA(512) 

INC=1 

OPEN ! UNI T-2 , NAME=NAME , TYPE= 'OLD ' , REC0RDSI2E= 1000 , 

X ACCESS = ' 1) I RECT ' , ASSOC I ATE VAR J ABLE= I NC ) 

READ!2'INC) (TDATA! I ) , 1=1 ,512) 

CLOSE ! UN I T =2 , D I SPOSE= ' SAVE ' ) 

RETURN 

END 


O o 

^ S3 

O ^ 
S pi 
«0 


SUBROUTINE PRONY ! ALMA , ALPH , ROOTR , ROOTI , W , N , M ) 

COMMON / PRO/ DELTA , DELT 
!.:OMMnN/NDATA/AVDATA ! 72 ) 

DIMENSION ALMA!N,N) ,ALPH!M),ROOTR(N) ,ROOTI(N) ,W(M) ,P0LEF!3S> 

COMPLEX RF.SMAT!3B,3R) ,POI E(3B) ,P0LTP(3B),RESIDU{36),CMP1 X,CI OB 
C 

SFT UP THE ALPHA MATRIX *#**»*»***»«**«*****iHHt»*«» 

f — * 

uj 

KD 


DO 10 1 = 1, N 


10 


K = I-1 

ALPH ( I ) =-AMDATA ( I+N ) 

DO 10 J=lrN 
ALtIA ( I , J ) =AUDATA ( J+K ) 

SOLUE FOR THE ALPHAS ^******tt*»^^n**it9^*** 

r. 

CALL S JM»( ALMA, ALPH rN^KS) 

IF{KS.£Q.O)GO TO 30 
TYPE B5 
GO TO 100 

30 ALPH(N+n = KO 

C 

SOLVE FOR POLES AS ZEROS OF Nth ORDER POLYNOMIAL ******* 

Q 

CALL PnLRT<ALPH,N,N»ROOTRrROOTJ rIER) % || 

IF( IER.EQ.O)GO TO 40 
TYPE B7,IER 
GO TO 100 
40 DO 22 1=1, N 

POLE < I ) =CMPLX ( ROOTR ( I ) , ROOTI ( I > ) 

22 POLE (I ) = ( 1 . /DELTA ) «CLOG C POLE ( J ) ) 

LOOK FOR POLES WITH POSITIVE REAL PARTS ♦»♦**»***** 

REMOVE WHEN FOUND **********************^ii* 

C 

l.=0 

DO 24 1=1, N 

IF{REAL<PQLE( I) ) .QT.O. )G0 TO 23 
L=L+J 
GO TO 24 

23 PQLE{I)=(0. ,0.) 

24 CONTINUE 

CK^f* COMPUTE RFSJDUES USING ONLY POLES WITH NEGATIVE REAL PARTS *** 


t-* 
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CALL flATRD ( RESMAT , POLE , RES IDU r POLEF » POLTP , L ) 

E5 FORMAT (JK. 'ALPHA MATRIX IS SINRULARl I ' ) 

67 FORMATnX, 'ERROR IN SUBROUTINE POLRT I ' r 5X » ' IER= ' , 15) 
100 RETURN 
END 


SUBROUTINE MATRD ( RESM » POLE r RESI , POLF r POLT , L > 

C 

C#^^*****^Hf^i«* ROUTINE TO SETUP MATRIX ANEI VECTOR FOR COMPUTATION OF 
C RESIDUES » AND TO PRINT OUTPUT »♦***♦**»*»*******»*♦ 

C0MM0N/DATA/RDATA(51Z) 

COMMON/ PARA/N , NN » NS , NBHI F 

COMMON/DAV/NAME,A 

COMMON/ PRO/DELTA » DELT 

DIMENSION RESOB) f POLF CL) rRCDATACSlZ) 

COMPLEX RESM(L,L) ,P0LE(36) ,RESI(L) • POLT(L) ,CEXPfCMPLX 
BYTE NAME(IG) ,A(1S) 

11 = 1 

DO 30 1=1, N 

IF(P0LE(I).EQ.(0.,0.))G0 TO 30 
POLTCII)=POLE(I) 

11 = 11+1 
30 CONTINUE 

DO JO 1=J.L 
K = I~1 

RES { I ) -RDATA ( 1 + { K^^NS ) ) 

DO 10 J=1,L 

1 0 RESM C I , J > =CEXP { POLT ( J ) *DELTA*K ) 

DO 40 1=J,L 

40 RFSl ( D=nMPLX(RES{I),0.0) 

C 

COMPUTE RESIDUES **««***«iHHt^Ht«****iHHt**** 

c 

CALL SIMQCCRESMfRESIrLfKS) 
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IF(KS.EQ.O)BO TO 70 
TYPE 69 

69 FORMATC DC, 'RESIDUE MATRIX IS SINGULAR!') 

GO TO 500 

C 

PRINT RESULTS ***#*»**»*****««*##**»♦*♦«♦»» 
C 

70 TYPE 100, NAME 
TYPE 110, A 
TYPE 400 

DO SO 1=1, L 

C IF(AIMAG{POLT(I) ) .LT.OIGO TO 80 

TYPE 4J0,PnLT(I) ,RESI(I) 

BO CONTINUE 
TYPE 412 

IF(DELT.EQ.50.E-12.QR.De.T.EQ.100.E-12)G0 TO 82 
SCAI-F = J ,B2399E-8 
FREQ=1 .591549E‘-7 
GO TO 84 

82 SCALE=9.70209E-10 

FRFG=8,4B5B8BE-9 
84 DO SO 1 = 1,1. 

IF(AIMAG(POLT(I)).LT.O)GO TO 90 
PQLFI I>=AIMAQ(POLT(I) }«FREQ 
POLTI I ) =POLT( I )*SCALE 
RESMAG=CABS(RESI ( I ) ) 

TYPE 411, POLT ( I ) , POLF (I ) , RESMAG 
90 CONTINUE 
C 

RECONSTRUCT NAMEFORM «**»♦*»***#**»* 

C 

C IF(NRFCnN,ER.O) GO TO 500 

DO 190 1=1,512 
RCDATA( I)=0.0 
DO 200 1=1, L 
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190 


lo 



200 PGLT( I)^P0LT( D/SCALE 
DO 250 1=1, L 

IF(AiriAG(POLT(I) ) .LT.OIGO TO 250 

IF(A)MAn(POLT( J ) ) .En,0)RESI (I)=C:i1PL>{(REAL(RESI{I) )».5, 

X AIMAG(RESKI))) 

IF(DELT.EQ.10.E-9)G0 TO 210 
NRFC=^00-NSHJF 
GO TO 230 
210 NREC=40 
230 DO 300 J=1,NREC 
•miF=CJ-3 )*I)ELT 

300 RCDATAf J)=RCDATA( J)+2«(EXP(REAL(P0LT{I) )*TIME) )* 

X (REALIRESJ ( J) J«COS{AIMAG{POI.T(I))*TIME)- 
X AltlAG(RESKI) )^tSIN( AIMAG(POLT( I ) )»TIME) ) 

250 CONTINUE 
C DO 310 1=1,50,5 

C3J0 TYPE 320,RCDATA(I ) ,RCDATA(I+1) ,RCDATAi( 1+2) ,RCDATA( 1+3) , 

C XRCDATACI+4) 

320 FORMAT ( IX, 5E 12, 4) 

C 

COMPUTE RMS ERROR -JHHHHHHHHHf* 

c 

RMSERR=0. 

RMSRSP=0. 

DO 350 I=1,NREC 

RMSERR =RMSERR+ ( RCDATA ( I ) -RDATA ( I ) ) **2 
350 RMSRHP=RM.SRSP+ (RDATA { I ) )**2 

RMSERR = { SQRT ( RMSERR/RMSRSP ) ) 3frlOO . 

TYPF 420 rNREC, RMSERR 
1 00 FORMAT { 26X , 1 GA 1 ) 

110 F0RMAT(27X, ISA) ) 

400 F0RMAT(20X, 'POLES', 35K, 'RESIDUES') 

410 FORMAT! 1X,E15.5,5X, '+J ' ,E15.5,5X,E15.5,5X, '+J',E15.5) 

411 FORMAT (5X.F) 0. 2, 5X, '+ J ' ,F1 0. 2 , 5X,FJ O. 2 , 10X,E) 5,5) 

412 F0RMAT<J5X, 'SCALED POLES', JOX, 'FRFRUENCVfMHZ ' , 9X, 'RESIDUE 
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420 F0RMAT(/lX, 'NREC=',I5,5X, 'RMS ERR0R=SEJZ,4, ' X') 
500 RETURN 
END 
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JSUBROimNF SI MR 
PLIRPnSE 

OBTAIN SOLUTION OF A SET OF SIMULTANEOUS LINEAR EQUATIONS, 
AX=B 


USAGE 

CALL SIMQ(A,B,N,KS) 

DESCRIPTION OF PARAMETERS 

A ~ MATRIX OF COEFFICIENTS STORED COLUMNWISE. THESE ARE 
DESTROYED IN THE COMPUTATION. THE SIZE OF MATRIX A IS 
N BY N. 

B - VECTOR OF ORIGINAL CONSTANTS (LENGTH N). THESE ARE 
REPLACED BY FINAL SOLUTION VALUES, VECTOR X. 

N - NUMBER OF EQUATIONS AND VARIABLES. N MUST BE .GT. ONE. 
KS - OUTPUT DIGIT 

G FOR A NORMAL SOLUTION 
1 FOR A SINGULAR SET OF EQUATIONS 


REMARKS 

MATRIX A MUST BE GENERAL. 

IF MATRIX IS SINGULAR , SOLUTION VALUES ARE MEANINGLESS. 

AN ALTERNATIVE SOLUTION HAY BE OBTAINED BY USING MATRIX 
INVERSION (MINV) AND MATRIX PRODUCT (GMPRDX. 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 
C 

C METHOD 

C METHOD OF SOLUTION IS BY ELIMINATION USING LARGEST PJVOTAI. 

C DIVISOR. EACH STAGE OF ELIMINATION CONSISTS DF INTERCHANGING 


^ f, . i /I e 
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C ROWS WHEN NECESSARY TO A^'QID DIVISION BY ZERO OR SMALL 

C ELRI1ENT5. 

C THE FORWARD SOLUTION TO OBTAIN VARJABI.F N IS DONE IN 

C N STAfiffS. THE BACK SOLUTION FOR THE OTHER VARIABLES IS 

C CALCULATED BY SUCCESSIVE SUBSTITUTIONS. FINAL SOLUTION 

C VALUES ARE DEVELOPED IN VECTOR B» WITH VARIABLE 1 IN B(l), 

C VARIABLE 2 IN B(2), VARIABLE N IN B(N>. 

C IF NO PIVOT CAN BE FOUND EXCEEDING A TOLERANCE OF O.Or 

C THE MATRIX IS CONSIDERED SINGULAR AND KS IS SET TO 1. THIS 

C TOLERANCE CAN BE MODIFIED BY REPLACING THE FIRST STATEMENT. 

C 

C 

c 

SUBROUTINE SIMQ ( A r B , N rKS ) 

DIMENSION A(1 )»B(1) 

FORWARD SOLUTION 

T0L=0.0 
KS-0 
JJ = -N 

DO 65 J=lrN 
JY=J+1 
JJ=JJ+N+1 
BIGA=0 
IT=JJ-J 
DO 30 I=J,N 

SEARCH FOR MAXIMUM COEFFICIENT IN COLUMN 
i: 

IJ=JT+I 

IF( ABS(BIGA)-ABS<A( IJ) ) ) 20,30,30 
20 BIGA=A(IJ) 

IMAX=I 
30 CONTINUE 
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TEST FOR PIVOT LESS THAN TOLERANCE (SINGULAR MATRIX) 

J F (ABS { B r GA ) -TOL ) 35 ,35. 40' 

35 KS=1 
RETURN 

INTERCHANGE ROMS IF NECESSARY 

40 I J =.1+N» (.1-2) 

IT=IMAX-J 
DO 50 K-JrN 
1 1 = I1+N 
I2=I1+IT 
SAVE=A( II ) 

A( II )=A{ I2> 

A( I2)=SAVE 

DIVIDE EQUATION BY LEADING COEFFICIENT 

50 A(IJ )=A( I J )/BIGA 
SAVE=B{IMAX) 

B(It1AX)^E( J) 

B( J)^SAVE/BIGA 

ELIMINATE NEXT VARIABLE 

IF(J-N) 55,70,55 
55 JQS^N*(.I-U 
DO 65 IX=JY,N 
IXJ=IQS+IX 
IT^J-IX 
DO BO JX=JYrN 
IXJX=N*(JX-1)+IX 
JJX=IXJX+IT 

•«4 
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EO A( D; JX)=A( IXJK)-(A(IXJ)*A( JJX) > 
G5 B( IK) ^B{ IX)-(B{J)*A( IXJ) ) 

BACK SOLUTION 

70 MY=N-J 
IT=N«N 

DO 80 J=lfNY 
IA=IT-J 
IB=N-J 
IC=N 

DO 80 K=1,J 

B( IB)=B{ IB)-A( IA)*B( IC) 

IA=IA-N 

80 rc=ic-i 

RETURN 

END 
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SUBROUTINE SJMRw 


i; 

C 

C 

C 


C 

C 

C 

C 

C 


C 

C 

C 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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PURPOSE 

OBTAIN SOLUTION OF A SET OF SIMULTANEOUS LINEAR EQUATIONS, 
AK=B, WITH COMPLEX COEFFICIENTS 


USAGE 

CALL SIMQC(A,B,N,KS) 


DESCRIPTION OF PARAMETERS ^ ,, 

A - MATRIX OF COEFFICIENTS STORED COLUMNWISE. THESE ARE 

DESTROYED IN THE COMPUTATION. THE SIZE OF MATRIX A IS "3 £ . 

N BY Hr g? 

B - VECTOR OF ORIGINAL CONSTANTS (LENGTH N), THESE ARE ^ 

REPLACED BY FINAL SOLUTION VALUES, VECTOR X. .© ^ 

N - NUMBER OF EQUATIONS AND VARIABLES. N MUST BE .GT. ONE. 

KS - OUTPUT DIGIT E 

0 FOR A NORMAL SOLUTION ^ 

1 FOR A SINGULAR SET OF EQUATIONS 


REMARKS 

MATRIX A MUST BE GENERAL. 

IF MATRIX IS SINGULAR , SOLUTION VALUES ARE MEANINGLESS. 

AN ALTERNATIVE SOLUTION MAY BE OBTAINED BY USING MATRIX 
INVERSION (MINV) AND MATRIX PRODUCT (GMPRD). 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 

METHOD 

METHOD OF SOLUTION IS BY ELIMINATION USING LARGEST PIVOTAL 
DIVISOR. EACH STAGE OF ELIMINATION CONSISTS OF INTERCHANGING 


V£) 


C ROWS WHEN NECESSARY TO AVOID DIVISION BY ZERO OR SMALL 

C ELEMENTS. 

C THE FORWARD SOLUTION TO OBTAIN VARIABLE N IS DONE IN 

C N STAGES. THE BACK SOLUTION FOR THE OTHER VARIABLES IS 

C CALCULATED BY SUCCESSIVE SUBSTITUTIONS. FINAL SOLUTION 

C VALUES ARE DEVELOPED IN VECTOR B, WITH VARIABLE 1 IN B(l), 

C VARIABLE 2 IN B(2)r VARIABLE N IN B(N). 

C IF NO PIVOT CAN BE FOUND EXCEEDING A TOLERANCE OF 0.0, 

C THE MATRIX IS CONSIDERED SINGULAR AND KS IS SET TO 1. THIS 

C TOl.ERANCE CAN BE MODIFIED BY REPLACING THE FIRST STATEMENT. 

C 

C 

c 

SUBROUTINE SJMflC( A,B,N,KS) 

COMPLEX A,B,BIGArSAVE 
DIMENSION A(l) ,B(1) 

C 

C FORNARD SOLUTION 

C 

T0L=0.0 
KS = 0 
JJ = -N 

DO G5 J=1 ,N 
JY=J+1 
JJ=JJ+N+1 
BIGA^^O 
IT=JJ-J 
DO 30 I=J,N 
C 

C SEARCH FDR MAXIMUM COEFFICIENT IM COLUMN 

C 

1J=IT+I 

IF(CABS(PJGA)-CABS«A( IJ) ) ) 20,30,30 
20 BIGA^A(IJ) 

IMAX=I 
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30 CONTINUE 

TEST FOR PIVOT LESS THAN TOLERANCE (SINGULAR HAT^^IX) 

IF (CABS(BIGA)-TOL) 35f35»40 
35 KS=1 
RETURN 

INTERCHANGE ROWS IF NECESSARY 


^0 lJ=J+N*(J-2) 

IT=IMAK“J 
DO 50 K^J.N 
Il=ri+N 
I2~ri+IT 
SAVE=A( II ) 

A(n)-A(I2) 

A(I2)=SAVE 

DIVIDE EQUATION BY LEADING COEFFICIENT 

50 A(I J )-A(IJ )/BIGA 
SAVF--B{ niAX) 

B{IMAX)-B( J) 

B( J)=SAVE/BIGA 

ELIMINATE NEXT VARIABLE 

IF(J-N) 55,70,55 
55 IOS==N*( ;-I) 

DO 65 IX=JY,N 
IXJ=IQS+IX 
IT=J-IX 
DO BO JX=JY,N 
IXJX=N»( JX-1 )+IX 
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JJX^IKJX+IT 

GO A(IXJK)=A(IXJ>I)-(A(IX*J)*A(JJX) ) 
(55 l?.( JX)=B(1X)-(B( J)«ACIKJ) ) 

C 

C BACK BflLllTION 

C 

70 NY=N-1 
IT=N*N 

DO 80 J^l.NY 
lA^IT-J 
IB=N-J 
IC=N 

DO BO K=1,J 

B(IB)=B( IB)-A(IA)*B<IC) 

IA=IA-N 
80 IC=IC-1 
RETURM 
END 
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c 

c 

i: 
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c 

c 

c 

n 

c 

c 

r. 

c 
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c 

c 

c 
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SUBROUTINE POLRT 
PURPOSE 

COMPUTES THE REAL AND COMPLEX ROOTS OF A REAL POLYNOMIAL 
USAGE 

- CALL POLRT(XCOFrCOF,M,ROOTR,ROOTI,IER) 

DESCRIPTION OF PARAMETERS 

XCOF -VECTOR OF M+1 COEFFICIENTS OF THE POLYNOMIAL 
ORDERED FROM SMALLEST TO LARGEST POWER 
COF -WORKING VECTOR OF LENGTH M+1 
M -ORDER OF POLYNOMIAL 

ROOTR-RESLILTANT VECTOR OF LENGTH M CONTAINING REAL ROOTS 
OF THE POLYNOMIAL 

ROOT r-RESULT ANT VECTOR OF LENGTH M CONTAINING THE 

CORRESPONDING IMAGINARY ROOTS OF THE POLYNOMIAL 
IFR -ERROR CODE NHFRE 
IER=0 NO ERROR 
IER=1 M LESS THAN ONE 
IER=2 M GREATER THAN 3S 

IER=3 UNABLE TO DETERMINE ROOT WITH 500 INTERATIONS 
ON 5 STARTING VALUES 
IER=4 HIGH ORDER COEFFICIENT IS ZERO 

REMARKS 

LIMITED TO 3BTH ORDER POLYNOMIAL OR LESS. 

FLOATING POINT OVERFLOW MAY OCCUR FOR HIGH ORDER 
POLYNOMIALS BUT WILL NOT AFFECT THE ACCURACY OF THE RESULTS. 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 


o c 

o 

O 

=■ 

-<ja 



c 
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c 

i: 

c 

i: 

c 

c 
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c 

c 

c 

c 

c 

c 

c 
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c 
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c 

c 

c 

c 

c 

c 


METHOD 


NEWTQN-RAPHSON ITERATIVE TECHNIQUE. 

ON EACH ROOT ARE PERFORMED USING THE 
RATHER THAN THE REDUCED POLYNOMIAL TO 
ERRORS IN THE REDUCED POLYNOMIAL. 


THE FINAL ITERATIONS 
ORIGINAL POLYNOMIAL 

avoid accumulated 


SHRROUT) NF POLRT ( XCOF , COF , M , ROOTR , ROOTI , I ER ) 
DIMENSION XCOFf J ) »COF( J ) r ROOTR ( 1 ) , ROOTI ( 1 ) 


IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS 

C IN COLUMN 1 SHOULD BE REMOVED FROM THF DOUBl E 
STATEMENT WHICH FOLLOWS. miooLfc 


DESIRED, THE 
PRECISION 


DOUBLE PRECISION XCOF, COF, ROOTR, ROOT I 


I HE C MUST ALSO BE REMOVED FROM DOUBLE PRECTSTOM QTflTPMirMTe 

rON'iTANT‘'lM'’qTfl^licl!!T'^^a®iS" MODIFIED BV CHflNOINB THE 

T STAIEMENT 78 TO l.OD-12 AND IN STATEMENT 12*^ TO 

cnsrip-EKFcin^liT^ precision results at^he 


IFIT-^0 

N=M 

IER-^0 




Ln 




n 

c 

c 


IF(XC0F(N+1) ) 10*25,10 
10 IF(N) 15,15,32 

SFT FRROR CODE TO J 

15 IFR=1 
20 RFTIIRN 
C 

C SET ERROR CODE TO 4 

C 

25 IER^4 
GO TO 20 
C 

C SET ERROR CODE TO 2 

C 

30 1FR=2 
GO TO 20 

32 IF(N-3B) 35,35,30 
35 NX=N 
NXX=N+1 
N2=l 

KJl = N+1 
DO 40 L-1,KJ1 
MT^KJl-L+1 
40 C0F(t1T)-XC0F(L) 

C 

C SFl INITIAL VAUJER 

C 

45 X0-. 005001 01 
Y0-- 0.01000101 
C 

C ZERO INITIAL VALUE COUNTER 

C 

1N = 0 
50 
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tn 

cn 



c 

C INCREMENT INITIAL VALUES AND COUNTER 

C 

Xn=-.10.*>»Y0 

Y0=-10.0*X 

C 

C SET X AND Y TO CURRENT VALUE 

C 

X=XO 
Y^YO 
IN=IN+1 
GO TO 50 
55 IFIT=1 
XPR=K 
YPR=Y 
C 

C EVALUATE POLYNOMIAL AND DERIVATIVES 

C 

50 JCT=0 
GO UX=0.0 
UY=0.0 
V =0.0 
YT=0.0 
XT=1.0 
U=C0F(N+1 ) 

IF<U) 65rl30,B5 
G5 DO 70 1=1, N 
L =N-I+1 
TEMP=COF(L) 

XT2=X^^XT-Y^^YT 

YTZ=X*YT+YitXT 

U=U+TEMP»XT2 

V=V+TEMP*YT2 

FI=I 

UX=UX+F Ii^XT*TEMP 
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UV-UY-FH^YT*TE«P 
>;T=XT2 
70 VT-=rT2 

SUMSQ-UX«UX+UY«UY 
IF(SUMSa) 75,110,75 
75 nx= ( o,”iJY-u«iix) /Hliiisa 
X=X+DX 

DY^- ( Uijy Y+V#UX ) /SUMSQ 
Y=Y+DY 

78 IF(DABS(DY)+DABS(DX)-1.0D-05) 100,^i0,80 

STEP ITERATION* COUNTER 

80 ICT::JCT+J 

IFdCT-500) 80,85,85 
85 IF{IFIT)100,a0,100 
90 IF{JN-5) 50,95,95 


SET ERROR CODE TO 3 


95 1ER=3 
GO TO 20 

100 DO 105 L=1,NXX 
HT-^KJl-L+1 
TEI1P=XCQF{MT) 
XCGF(MT}=COF(L) 

105 CQF(L)=TEHP 
ITEMP-N 
N=NX 

NX=ITEMP 

IF(IFIT) 120,55,120 
110 IF(IFIT) 135,50,115 

nr; x=xpr 

Y = YPI? 

120 IFIT=0 
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122 IF(DABS(Y)-1 .0D-4*DA8S(X> ) 135,125,125 
125 AI.PHA=X+X 

SUMSa=X«K+Y»Y 
H^N-2 
GO TG 140 
130 >C=0.0 
NK=NX-1 
NXX=NXK-1 
135 Y=0.0 

SUMSQ^O.O 
ALPHA =X 
N-N-1 

140 CQF{2)=C0F(2>+ALPHA#C0F( 1) 

145 DO 150 L=2,N 

150 i:OF{|_+l )=CDF(L+l)-f-ALPHA«(:OF(L)-SUMSG«COF(L-i) 
155 RD0TI(N2)=Y 
R00TR{N2)=K 
H2-H2+1 

IF I SUHSQ ) ISO , 1B5 , ISO 

mo 

GO TO 155 
165 IF4NJ 20,20,45 
Bm 
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